Gene-regulatory enhancers have been identified using various approaches, including evolutionary conservation, regulatory protein binding, chromatin modifications, and DNA sequence motifs. To integrate these different approaches, we developed EnhancerFinder, a two-step method for distinguishing developmental enhancers from the genomic background and then predicting their tissue specificity. EnhancerFinder uses a multiple kernel learning approach to integrate DNA sequence motifs, evolutionary patterns, and diverse functional genomics datasets from a variety of cell types. In contrast with prediction approaches that define enhancers based on histone marks or p300 sites from a single cell line, we trained EnhancerFinder on hundreds of experimentally verified human developmental enhancers from the VISTA Enhancer Browser. We comprehensively evaluated EnhancerFinder using cross validation and found that our integrative method improves the identification of enhancers over approaches that consider a single type of data, such as sequence motifs, evolutionary conservation, or the binding of enhancer-associated proteins. We find that VISTA enhancers active in embryonic heart are easier to identify than enhancers active in several other embryonic tissues, likely due to their uniquely high GC content. We applied EnhancerFinder to the entire human genome and predicted 84,301 developmental enhancers and their tissue specificity. These predictions provide specific functional annotations for large amounts of human non-coding DNA, and are significantly enriched near genes with annotated roles in their predicted tissues and lead SNPs from genome-wide association studies. We demonstrate the utility of EnhancerFinder predictions through in vivo validation of novel embryonic gene regulatory enhancers from three developmental transcription factor loci. Our genome-wide developmental enhancer predictions are freely available as a UCSC Genome Browser track, which we hope will enable researchers to further investigate questions in developmental biology.
Eukaryotic gene expression is regulated by a highly orchestrated network of events, including the binding of regulatory proteins to DNA, chemical modifications to DNA and nucleosomes, recruitment of the transcriptional machinery, splicing, and post-transcriptional modifications. Enhancers are genomic regions that influence the timing, amplitude, and tissue specificity of gene expression through the binding of transcription factors and co-factors that increase transcription (as reviewed in [1], [2]). In humans, genetic variation in enhancer regions is implicated in a wide variety of developmental disorders, diseases, and adverse responses to treatments [3], [4], [5]. Enhancers have been discovered in introns, exons, intergenic regions megabases away from their target genes [6], and even on different chromosomes [7]. An enhancer frequently drives only one of many domains of a gene' s expression [8], [9] and different cell types accordingly exhibit considerable differences in their active enhancers [10], [11]. This modularity enables the creation of complex regulatory programs that can evolve relatively easily between closely related species [12], [13]. Individual enhancers were initially identified using transgenic assays in cultured cell lines [14], [15] and later in vivo in model organisms, such as mouse, Drosophila, and zebrafish. In the in vivo experiments, a construct containing the sequence to be tested for enhancer activity, a minimal promoter, and a reporter gene (e. g. , lacZ) is injected into fertilized eggs, and transgenic individuals are assayed for reporter gene expression. Early efforts to find enhancers at the genome scale used comparative genomics. Several studies assayed non-coding regions conserved across diverse species for enhancer activity [16], [17], [18], since functional non-coding regions likely evolve under negative selection. This approach identified many enhancers at a range of levels of evolutionary conservation [19], [20], [21]. However, relying on evolutionary conservation alone has several shortcomings: many characterized enhancers are not conserved between species [22], non-coding conservation is not specific to enhancer elements, and evolutionary patterns provide little information about the tissue and timing of enhancer activity. Enhancer prediction has been revolutionized by recent technological advances, including chromatin immunoprecipitation coupled with high-throughput sequencing (ChIP-seq) [23], RNA sequencing (RNA-seq), and sequencing of DNaseI-digested chromatin (DNase-seq) [24] or formaldehyde-assisted isolation of regulatory elements (FAIRE-seq) [25]. These “functional genomics” assays enable genome-wide measurement of histone modifications, binding sites of regulatory proteins, transcription levels, and the structural conformation of DNA. The ENCODE project [26], FANTOM project [27], and similar studies focused on specific cell types [28], [29] have dramatically increased the amount of publicly available functional genomics data. Functional genomics studies revealed several genomic signatures of active enhancers. For example, known enhancers are associated with the unstable histone variants H3. 3 and H2A. Z [30], [31] and low nucleosome occupancy [32], although these chromatin states are not unique to enhancers. Monomethylation of lysine 4 on histone H3 (H3K4me1), a lack of trimethylation at the same site (H3K4me3), and acetylation of lysine 27 on histone H3 (H3K27ac) may distinguish active enhancers from promoters [10], [33], [34], enhancers that are “poised” for activity later in development [35], [36], and regulatory elements that repress gene expression [37], [38]. Additional features that pinpoint specific classes of active enhancers include binding of the transcriptional cofactor p300/CBP [18], [39], [40], [41], clusters of transcription factor (TF) binding sites [42], [43], [44], [45], and enhancer RNA transcription (eRNAs) [46]. Collectively, functional genomics data have pinpointed the locations of many novel enhancers and yielded insights into sequence and structural determinants of enhancer activity. However, these patterns have not proven to be universal [47], [48], and there is unlikely to be a single chromatin signature that identifies all classes of enhancers [11], [49], [50]. Given the complexity of these functional genomics data sets, computational methods have been developed to improve and generalize the enhancer predictions made from simple combinations of these data. Support vector machines (SVMs) and linear regression models trained to interpret DNA sequence motifs underlying known enhancers have successfully identified novel enhancers active in heart [51], hindbrain [52], and muscle [53] development. Another approach used SVMs to learn patterns of short DNA sequence motifs that distinguish markers of potential enhancers, such as p300 and H3K4me1, in different cellular contexts [54], [55]. Random forests have been used to predict p300 binding sites from histone modifications in human embryonic stem cells and lung fibroblasts [56]. Machine-learning algorithms have also been applied to the related problem of selecting functional TF binding sites out of the thousands of hits to a TF' s binding motif throughout the genome [57], [58], [59], [60], [61], [62], [63]. Finally, two groups have taken a less supervised approach and used hidden Markov models (ChromHMM) [64] and dynamic Bayesian networks (Segway) [65] to segment the human genome into regions with unique signatures in ENCODE data and then assigned potential functions, such as enhancer activity, to these states. While rich datasets coupled with sophisticated algorithms have successfully identified many novel enhancers, comprehensive enhancer prediction is challenging for two main reasons. First, no single type of data is currently sufficient to identify all enhancers active in a given context. Many of the approaches described above use a single mark or motif as a proxy for an enhancer, but this gives an incomplete representation of all biologically active enhancers. Second, while a great deal of functional genomics data are available for different cell lines and tissues, it is not understood how informative experiments in a given cellular context are indicative of enhancer activity in other contexts. With these issues in mind, we introduce EnhancerFinder, a new two-step machine-learning method for predicting enhancers and their tissue specificity. In machine learning, a classification algorithm is trained to distinguish between labeled training examples (e. g. , enhancers and non-enhancers) based on features of these labeled examples (e. g. , evolutionary conservation, chromatin signature, DNA sequence). The trained classifier can then be used to predict the labels for uncharacterized genomic regions (e. g. , which ones are enhancers). Our approach employs two rounds of a supervised machine-learning technique called multiple kernel learning (MKL) [66], [67]. MKL is based on the theory of SVMs [68], but provides greater flexibility to combine diverse data (e. g. , evolutionary conservation, sequence motifs, and functional genomics data from different cellular contexts) and to interpret their relative contributions to the resulting predictions. Our implementation of EnhancerFinder applies MKL in two steps with the goal of generating a genome-wide set of developmental enhancers to better characterize gene regulation during development. The algorithm, which is trained using in vivo validated enhancers from the VISTA enhancer database [69] and publicly available genomic data, first aims to distinguish human developmental enhancers from the genomic background and then in a second step predicts enhancer tissue specificity. In contrast to most other enhancer prediction strategies, which are trained on epigenetic marks or sequence motifs that serve as a proxy for a subset of all active enhancers, our use of a heterogeneous and in vivo validated set of enhancers, enables us to investigate the complex suite of features that underlie active regulatory regions. With appropriate training data, EnhancerFinder could be applied to study gene regulation at other developmental stages. Our analyses demonstrate that EnhancerFinder' s integration of diverse types of data from different cellular contexts significantly improves prediction of validated enhancers over approaches based on a single context or type of data. We find that enhancers active in some developmental contexts are easier to identify than others. Applying EnhancerFinder to the entire human genome allowed us to predict more than 80,000 developmental enhancers, with tissue-specific predictions for brain, limb, and heart. These predictions significantly overlap known non-coding regulatory regions and are enriched near relevant genome-wide association study (GWAS) lead single nucleotide polymorphisms (SNPs) and genes expressed in the predicted tissue. To illustrate the utility and accuracy of our genome-wide enhancer predictions, we used them to investigate the enhancer landscape near three developmentally expressed genes. First, we screened predicted enhancers near FOXC1 and FOXC2 in transgenic zebrafish, and found that 70% (7 of 10) of tested EnhancerFinder predictions have confirmed (6) or suggestive (1) developmental enhancer activity. In addition, we validated a novel cranial nerve enhancer near the ZEB2 locus using a transgenic mouse enhancer assay. Taken together, our results suggest that the EnhancerFinder approach of integrating diverse data sets significantly improves prediction of biologically active enhancers, providing high-confidence candidate enhancers for studies in developmental gene regulation.
Step 1 of our pipeline aims to distinguish all enhancers active in the context of interest (e. g. , a specific developmental stage) from non-enhancer regions. Step 2 then builds classifiers to predict the tissues in which the enhancer candidates from Step 1 are active. This two-step approach allows us to accurately identify enhancers, while also distinguishing their tissues of activity. We train and evaluate EnhancerFinder using the VISTA Enhancer Browser, which at the time of our analysis contained over 700 human sequences with experimentally validated enhancer activity in at least one tissue at embryonic day 11. 5 (E11. 5) in transgenic mouse embryos. VISTA also contained a similar number of regions without enhancer activity in this context. E11. 5 in mouse development roughly corresponds to E41 (Carnegie stage 17 [70]) in human development. In Step 1 of EnhancerFinder, we used all 711 VISTA enhancers as positive training data, and for negative training data, we created a set of 711 random regions matched to the length and chromosome distribution of the positives to represent the genomic background. We did not use the VISTA negatives as negative training examples in Step 1, because they are not representative of all non-enhancer regions (see below). Our goal in Step 1 is to develop a method that can be used to scan the whole genome and distinguish developmental enhancer regions from non-enhancer regions. The second step of EnhancerFinder aims to distinguish enhancers active in a given embryonic tissue from non-enhancers and enhancers active in other tissues. We consider all enhancers in VISTA with activity in a tissue of interest as positives and all other regions in VISTA (including regions not active at E11. 5) as negatives (see Methods). This second step that includes enhancers active in other tissues as negatives in the training proves to be essential for obtaining high specificity in predicting tissue of activity (see below), and it is important to do this in two steps rather attempting to distinguish enhancers of a given tissue from genomic background in one step. To evaluate EnhancerFinder, we compared it to several commonly used enhancer prediction approaches. Unless otherwise noted, we evaluated the performance of all prediction algorithms using 10-fold cross validation to compute the area under the curve (AUC) for receiver operating characteristic (ROC) curves. We also computed precision-recall curves (Figure S1) and compared power at a low false positive rate. Because EnhancerFinder learns enhancer signatures from a training data set, we first explored biases in the VISTA enhancers that might affect how well EnhancerFinder could generalize to the whole genome. The genomic regions tested by VISTA were not selected randomly, and thus their positives do not represent a random sample of active enhancers. Nearly all regions tested by VISTA are evolutionarily conserved across mammals (706 of 711 positives and 727 of 736 negatives). Since our goal is to predict a broadly applicable, high confidence set of developmental enhancers, we did not include this feature when making genome wide predictions. However, with this bias in mind, we did evaluate several models that incorporate the degree of evolutionary conservation (see below). In addition to conservation, several studies deposited in VISTA have considered enhancer-associated proteins and histone marks, such as p300, H3K27ac, and H3K4me1. We collected all data sets of these types from ENCODE and computed their overlap with VISTA enhancers. Fewer than half of the VISTA positives are marked by all three of p300, H3K27ac, and H3K4me1 (from any data set), with substantial percentages marked by only one or two and 13% (93/711) marked by none (Figure S2). These findings indicate that VISTA positives are not highly biased towards a single type of ChIP-seq feature, motivating us to include these features in our genome-wide predictions, with the caveat that the trends we observe for VISTA positives might not generalize to all classes of enhancers. Our analysis also suggests that the standard practice of equating active enhancers with all regions marked by a single ChIP-seq feature, or even the union of overlapping peaks from several ChIP-seq experiments, will fail to identify all active enhancers in a given context. EnhancerFinder predicts enhancers by integrating classifiers based on distinct data types. In our first evaluation of EnhancerFinder, we consider: functional genomics data, evolutionary conservation patterns, and DNA sequence motifs. Combining these different approaches enables EnhancerFinder to accurately distinguish enhancers from the genomic background (Figure 2A; AUC = 0. 96). The functional genomics component of EnhancerFinder (which we refer to as All Functional Genomics) is a linear SVM that incorporates 2469 datasets generated by the ENCODE project and smaller scale studies. These include DNaseI hypersensitivity data and ChIP-Seq for p300, many histone modifications, and many TFs from many adult and embryonic tissues and cell lines (Table S1). DNA sequence patterns are integrated via a 4-spectrum kernel (DNA Motifs), which summarizes the occurrence of all length four DNA sequences (4-mers) in input regions [71]. We found that little was gained by increasing k, considering multiple k simultaneously, or incorporating knowledge of transcription factor binding site (TFBS) motifs as in a previous approach [52] (Figures S3 and S4). Finally, evolutionary conservation information is incorporated with a linear SVM that uses mammalian phastCons scores [72] as features (Evolutionary Conservation). One motivation for developing EnhancerFinder was to explore whether combining previous successful approaches to enhancer prediction would improve performance. Each of the classifiers combined in EnhancerFinder is representative of a different strategy for predicting enhancers. Thus, we compared the performance of EnhancerFinder to each of its constituents, which are SVMs trained on the same enhancer data as EnhancerFinder, but using only one type of the data features (e. g. , only sequence motifs). EnhancerFinder significantly outperformed each of the individual classifiers (Figure 2A; p = 2. 0E-7 for Evolutionary Conservation, p = 2. 6E-8 for DNA Motifs, and p = 4. 4E-16 for All Functional Genomics, McNemar' s test), suggesting that these different types of data capture unique aspects of enhancers that are not completely encompassed by any single data type. Not surprisingly, we found that of the three component classifiers in EnhancerFinder, Evolutionary Conservation yields the best performance (AUC = 0. 93). As noted above, nearly all regions tested for enhancer activity by VISTA (positives and negatives) are evolutionarily conserved compared to the genomic background. Nonetheless, considering additional features significantly improved predictions. The DNA Motifs (AUC = 0. 88) and All Functional Genomics (AUC = 0. 89) classifiers also exhibit strong performance, but also do not perform as well as the combined classifier. EnhancerFinder has nearly twice the power of any of the individual classifiers at a 5% false positive rate (FPR), and its power advantage is even larger at lower FPRs. All Functional Genomics, DNA Motifs, and Evolutionary Conservation achieve roughly similar performance from different feature data, but each individual classifier predicts a somewhat different set of enhancers during evaluation (Figure 2B). Roughly two-thirds of the enhancer predictions are shared between the three classifiers. The improvement provided by combining these data argues that these data sources are indeed complementary. We also compared EnhancerFinder' s performance with several current computational methods used to identify enhancers. We were able to make the most direct comparison with CLARE, a popular method for identifying enhancers from DNA sequence data, i. e. , transcription factor binding site motifs and other sequence patterns [73]. This approach, which has been successfully applied in several contexts [51], [52], [53], [74], makes few assumptions about the input, and is publicly available as a web server. On our Step 1 enhancer prediction task, we find that CLARE achieves an ROC AUC of 0. 79. This is much lower than DNA Motifs (AUC = 0. 88), our approach based on sequence data alone, and the full EnhancerFinder (AUC = 0. 96; Figure 2C). At a 5% FPR, the power of CLARE is about 20%, compared to approximately 30% for DNA Motifs and more than 60% for EnhancerFinder. Comparisons with additional methods were complicated by the fact that most were developed in different contexts. We designed EnhancerFinder specifically to predict biologically active developmental enhancers. Most existing approaches focus on data from a single cell line and define enhancers based on specific enhancer-associated marks or proteins (such as p300 in human embryonic stem cells) rather than biological activity. Thus, we did not anticipate that they would perform as well as EnhancerFinder at developmental enhancer prediction. However, since the predictions of these methods are commonly used outside the specific contexts in which they were made, we believe that it is useful to evaluate how well they can identify developmental enhancers and how much the EnhancerFinder approach applied to developmental enhancers improves on their performance. In particular, we compared EnhancerFinder to ChromHMM and Segway [64], , two unsupervised machine learning methods for segmenting the genome into a small number of functional “states” based on consistent patterns in ENCODE data for individual cell lines. The states resulting from the segmentations of each cell line' s data are annotated by hand into predicted functional classes, which include enhancer activity. To evaluate these methods, we considered the states overlapping our training and testing regions. Any region with an overlapping enhancer state was considered a predicted enhancer and all others were predicted non-enhancers. In this way, we obtained a single point in ROC space for the state predictions. Since there is no score or confidence value associated with the state assignments, a full ROC curve could not be created for these methods. Figure 2C gives the performance for several versions of ChromHMM and Segway based on ENCODE data from different cell lines. Both methods perform better than random, but considerably worse than EnhancerFinder and CLARE (p≈0). We stress that, in contrast to our supervised method, these methods were not explicitly trained to perform the same task as EnhancerFinder, and thus we did not expect them to perform as well as EnhancerFinder. Indeed, these results argue that their utility in identifying developmental enhancers is limited compared to specialized approaches. As illustrated above, our machine learning prediction and evaluation framework enabled us to quantitatively explore the utility of different genomics datasets in enhancer prediction by creating classifiers based on different types of data (i. e. , sequence motifs, evolutionary conservation, and functional genomics) and comparing their performance. We also used this framework to investigate other questions about the utility of different subsets of these data for enhancer prediction. For example, one might expect that some of the datasets included in All Functional Genomics (e. g. , experiments in cancer cell lines or adult tissues) would not be as useful as others (e. g. , experiments in embryonic tissues) for predicting developmental enhancers, and that limiting the features examined by the classifier to the most relevant experiments might improve performance. To explore this hypothesis, we trained linear SVM classifiers to predict VISTA enhancers (as in Step 1 of EnhancerFinder) based on different subsets of all the functional genomics features (Table 1) and compared their performance. First, we considered a collection of 244 datasets from embryonic tissues and cell lines (Embryonic Functional Genomics). Next, we created a classifier that considers data from a wider range of contexts by training a linear SVM using a large, manually curated set of 509 potentially relevant functional genomics data sets (Relevant Functional Genomics). This set includes embryonic datasets, along with additional DNaseI and ChIP-seq data from adult tissues and cell lines related to the dominant tissues of activity in VISTA. For example, we included data from human cardiac myocytes, since there are many developmental heart enhancers in our training examples. We compared these to the All Functional Genomics classifier described above that uses all 2496 functional genomics features. All Functional Genomics (AUC = 0. 89) performed slightly, but not significantly, better than Relevant Functional Genomics (AUC = 0. 87; p = 0. 16), and both significantly outperformed Embryonic Functional Genomics (AUC = 0. 83; p = 9. 2E-9 and p = 2. 7E-6, respectively) (Figure 3A). At low FPRs, the differences in power between these classifiers were modest. The Embryonic Functional Genomics classifier included the most time-appropriate datasets, yet its performance was improved by including additional data sets that seem less relevant to our classification problem a priori. Thus, we conclude that it can be advantageous to consider a range of functional genomics features, especially when few features are available from the context of interest. The utility of these additional data sets might indicate that some enhancer features are stable across cell types and developmental stages, but it could also reflect information these data provide about genomic regions that are not active enhancers during development (see Discussion). We also explored the utility of individual functional genomics datasets that are often used as proxies for developmental enhancers by creating three linear SVM classifiers: H3K27ac, H3K4me1, and p300. These SVMs were trained to distinguish VISTA positives from the genomic background (Step 1) using all available data of the specified type from ENCODE, which include a range of cell types and tissues (Table S1). All three classifiers performed better than random (Figure 3B). H3K4me1 (AUC = 0. 72) and p300 (AUC = 0. 68) performed similarly (p = 0. 25), with p300 performing best at low FPRs and H3K4me1 best at higher FPRs. Both significantly outperformed H3K27ac (AUC = 0. 61; p = 9. 4E-15 and p = 5. 5E-9, respectively); however, we caution against extrapolating from this comparison, since it may reflect biases in the feature sets available and the VISTA positives. Since combinations of these features are often used to predict enhancers, we next trained a linear SVM classifier (Basic Functional Genomics) that includes all three data types together. The combined classifier significantly outperforms all the individual classifiers (AUC = 0. 77; p<2E-7 for each), suggesting that each data type contributes unique information about enhancer activity. Also, all four SVM classifiers achieved much better performance than the common approach of simply considering regions overlapping with these data (Figure S5). EnhancerFinder also learns weights for individual features within classifiers that reflect their contribution to the enhancer predictions. We found that features known to be associated with enhancer activity in relevant cellular contexts generally receive positive weights, while those associated with other types of elements received negative weights (Text S1 and Figure S6). In the previous sections, we focused on generic developmental enhancer prediction (Step 1 of EnhancerFinder). Step 2 of EnhancerFinder applies a second round of MKL to refine and further annotate predicted enhancers from Step 1 (Figure 1). In this study, Step 2 consists of training an MKL classifier to distinguish VISTA enhancers active in a given tissue from VISTA regions without activity in that tissue, i. e. , non-enhancers from VISTA plus enhancers for other tissues. We did not require that the positive training examples be active only in the tissue of interest. Using the same feature data as in Step 1, we created tissue-specific classifiers for all tissues with more than 50 examples in VISTA: forebrain, midbrain, hindbrain, heart, limb, and neural tube. The performance of EnhancerFinder' s tissue specificity predictions varied dramatically between tissues (Figure 4), with the best performance for heart (AUC = 0. 85), followed by limb (AUC = 0. 74), forebrain (AUC = 0. 72), midbrain (AUC = 0. 72), hindbrain (AUC = 0. 69), and neural tube (AUC = 0. 62), which was the worst of the tested tissue classifiers, but better than random. We combined all brain enhancers into one class, and the performance of this generic brain classifier was similar to that of the more specific brain classifiers (AUC = 0. 73). The EnhancerFinder tissue-specific classifiers trained with all data types performed well for most tissues (Table 1); however, classifiers based on functional genomics alone often performed as well as the full EnhancerFinder classifier, suggesting functional genomics data are more informative about developmental enhancer tissue specificity than degree of conservation or sequence motifs. Most previous efforts to predict tissue-specific enhancers have performed a single training step using enhancers or enhancer marks present in the tissue of interest as positives and non-enhancer regions or the genomic background as negatives. To test whether our two-step method improves upon these previous approaches, we trained one-step MKL tissue-specific classifiers and compared their predicted tissue distributions to those of validated enhancers from the VISTA database (Figure 5A). First, we trained a set of tissue-specific classifiers using enhancers active in each tissue as positives and the genomic background as negatives. These classifiers predict very similar sets of enhancers regardless of the target tissue; and they vastly overestimate the number of enhancers that are active in multiple tissues (95% of predictions versus 8% of VISTA) and the number of true enhancers of each tissue (Figure 5B). In contrast, classifiers trained as in Step 2 of EnhancerFinder, i. e. , using tissue-specific enhancers as positives and a mix of enhancers active in other tissues and regions with no activity in VISTA as negatives, show much greater tissue-specificity in their predictions (76%) and a similar amount of overlap as among known enhancers (Figure 5C). The relative ease of identifying heart enhancers is likely due to several unique characteristics. Known heart enhancers at E11. 5 are more evolutionarily conserved than genomic background, but significantly less conserved than enhancers in other tissues [39], [41]. In addition, we observed that heart enhancers at this developmental stage are uniquely close to the nearest transcription start site (TSS) (Figure S7). These two patterns are consistent with a recent study of mouse enhancers from different developmental stages [75]. Finally, we observed that E11. 5 heart enhancers have an unusually high GC content (49%) compared to enhancers of other tissues at E11. 5 (∼40%). A simple classifier based solely on the GC content of a region performs nearly as well as our full classifier for heart enhancers (Figure S8). In contrast, sequence-based classifiers do not perform well on the other tissues whose enhancer GC content is not significantly different from the genomic background (Table 1). The high GC content of heart enhancers is not due to overlap with CpG islands. Only about 4% of VISTA enhancers overlap with a CpG island, and this number is consistent across tissues. We also did not find enrichment for any known GC-rich transcription factor binding site motifs in VISTA heart enhancers. We do see, however, that repeat regions in heart enhancers are depleted for the very AT-rich repeats seen in other enhancers, and that most of the repeat regions in heart enhancers are 40–60% GC. Our results suggest the possible existence of unknown GC-rich motifs that may be important for gene regulation in the cardiac lineage. The heart classifier based on functional genomics data alone exhibits strong performance compared to other tissue-specific classifiers as well (Table 1). It is possible that this is due to the presence of feature data from contexts more relevant to developmental heart activity than to other tissues, rather than unique attributes of the heart enhancers themselves. Indeed, the highest weighted features in the heart functional genomics classifier come from heart tissues. However, the performance of the heart classifier based only on functional genomics data does not decrease substantially when we exclude data from the most relevant contexts: embryonic heart tissue, adult hearts, and stages of a directed differentiation of stem cells into cardiomyocytes (ROC AUC = 0. 85). Thus, it is possible that feature data from less obviously relevant contexts are more informative about heart activity than for other tissues. We suspect that the ease of distinguishing heart enhancers may be due to the earlier development of the heart compared to other tissues (see Discussion). One of the main motivations for developing algorithms that can distinguish active enhancers is to apply them to unannotated genomic regions to aid the exploration and interpretation of the gene regulatory landscape of the human genome (Figure 1). To produce a genome-wide set of candidate developmental enhancers, we divided the genome into 1. 5 kb blocks overlapping one another by 500 bp and applied Step 1 of EnhancerFinder to each of these regions. EnhancerFinder produces a score for each region; positive scores indicate membership in the positive set (enhancers), and negative scores indicate membership in the negative set (non-enhancers). To focus on high confidence predictions in this genome-wide analysis, we used the cross-validation-based evaluation described above to find a 5% FPR score threshold, and only considered regions exceeding this threshold. After merging overlapping positive predictions, we identified 84,301 developmental enhancers across the human genome with median length of 1,500 bp and total genome coverage of 183,695,500 bp (5. 86%). The 5% FPR threshold we used corresponds to a 65% true positive rate (TPR). To calculate the false discovery rate (FDR), we must estimate the unknown fraction of 1. 5 kb blocks of the human genome that harbor developmental enhancer regions. If this fraction were as high as 50%, a 5% FPR would correspond to a 9% FDR. If instead we estimate that 10% of 1. 5 kb windows contain a developmental enhancer, we see an FDR of 47% at a 5% FPR. While this may seem high, our recent analysis of predicted enhancers with human-specific substitution rate acceleration found a lower failure rate at E11. 5 (17%, 5/29) [74], and only three of ten tested predictions did not validate with confirmed or suggestive activity in our zebrafish assay (see below). This suggests that the FDR may be lower in experimental applications, especially when predicted enhancer regions are analyzed in the context of other relevant data. However, to accurately measure the true FDR would require experimental testing of a very large, random set of EnhancerFinder predictions, which is beyond the scope of this study. In our genome-wide analysis, we used the smaller Relevant Functional Genomics data set in order to reduce the computational time required. We also did not include evolutionary conservation data, because the positives in our training data are almost universally conserved. While most enhancers likely exhibit some evolutionary conservation, this extremely high fraction is likely due to bias in the selection of the tested regions in VISTA and could reduce our ability to detect less highly conserved novel enhancers genome-wide (see Discussion). The resulting conservation-free classifier still performed extremely well in cross validation (AUC = 0. 92). Supporting this approach, non-conserved regions make up over 20% of our genome-wide enhancer predictions. As noted above, we did not observe any other dramatic biases in the feature data associated with human VISTA enhancers. Next, we applied Step 2 of EnhancerFinder to all enhancer regions predicted in Step 1. We focused on brain, limb, and heart, because these tissues are highly represented in VISTA and have been extensively studied in previous analyses of developmental enhancers. We predicted 7,400 limb enhancers, 19,051 heart enhancers, and 11,693 brain enhancers (Figure 6) at a 5% FPR threshold tuned separately for each tissue. Since EnhancerFinder makes predictions for each tissue independently, there are no constraints on the distribution of tissues in the resulting genome-wide predictions. Nonetheless, we find a high level of tissue-specificity; nearly 90% of the limb, heart, and brain enhancers are predicted to be active in just one of the three tissues. All genome-wide enhancer predictions are available as tracks for import into the UCSC Genome Browser (Data File S1). These lists of high-confidence tissue-specific enhancers should not be viewed as exhaustive; we found thousands of regions with positive, but less significant scores from Step 2 of EnhancerFinder. To characterize and further validate our genome-wide enhancer predictions, we examined their genomic distribution with respect to several independent indicators of function (details in Text S1). Genes near brain and heart enhancers are enriched for expression in relevant tissues (Tables S2 and S3). Similarly, Gene Ontology (GO) Biological Process enrichment analyses of nearby genes suggest that our predicted developmental enhancers target genes that function in relevant cell types and tissues (Figure 6). The most prevalent transcription factor binding site motifs found in the sequences of predicted enhancers differed between enhancers of different tissues and included many relevant developmental TFs (Table S4). Finally, our predicted enhancers contain 676 lead SNPs associated with significant effects in GWAS (Table S5); this is significantly more than expected at random (permutation p<0. 001). Taken together, these analyses suggest that EnhancerFinder identifies many active regulatory regions that contain functionally relevant variation. Our tissue-specific enhancer predictions give valuable annotations to thousands of non-coding regions of the human genome that had not previously been linked to developmental regulation. For example, thousands of SNPs associated with disease by GWAS are in non-coding regions with limited functional annotations [76]. Our genome-wide enhancer predictions provide a resource for exploring the mechanisms and functional effects of these uncharacterized GWAS hits. To demonstrate that genome-wide EnhancerFinder predictions can facilitate the discovery of functional regulatory elements, we present two case studies in which we identify and validate novel enhancers near genes active during development.
We chose to define developmental enhancers for training as genomic regions that are experimentally shown to activate gene expression in vivo in embryonic mouse assays. We believe that this definition is better suited to identifying regions for further exploration and experimental characterization than approaches based on single data sources, such as p300, H3K4me1, or H3K27ac, associated with enhancers in individual cell lines. We showed that our predicted enhancers, based on this biologically active definition, significantly overlap data sets commonly used as proxies for enhancer activity, such as H3K27ac and p300 binding. However, these other data alone are not sufficient to identify all enhancers, as we demonstrated for H3K27ac, H3K4me1, and p300 in Figure 3B. Similarly, when we evaluated the ability of other computational methods to identify enhancers, we find that they perform better than random, but that EnhancerFinder significantly outperforms them at identifying biologically active developmental enhancers. This is not surprising given the different contexts in which some enhancer predictions, such as those from ChromHMM and Segway, were developed. While EnhancerFinder could be used to predict enhancers in well-characterized cell lines, it is particularly useful at identifying enhancers in complex tissues that contain multiple cell types and in cell types that do not have much specific functional genomics data available. Other computational approaches to enhancer prediction have focused on identifying enhancers in individual cell types using functional genomics data from the same cells [56] or using the differences in cell type specific transcription factor binding to identify cell-type specific binding motifs [61]. These methods generally perform well, but they do not address enhancer prediction in cell types with little or no functional genomics data, or in tissues that contain multiple cell types. Data such as p300 binding sites and H3K4me1 have been used in previous studies to identify enhancers, and these data are major contributors to our enhancer predictions. However, data from other sources and contexts less directly associated with enhancer activity provide complementary information that improves our predictions. Some of these data may be negatively correlated with enhancer activity, allowing EnhancerFinder to learn what features distinguish regions that are not developmental enhancers. Others may help reinforce patterns present in data from more relevant contexts, reflecting some degree of stability in the features of enhancer regions across developmental stages and cell types. For example, we found that features measured in embryonic stem cells are quite useful for E11. 5 enhancer prediction; their removal from the classifier degrades performance and/or they have large (positive or negative) MKL weights. Examination of these features suggests that some identify “poised” regions that will become active enhancers upon differentiation, while others seem to help distinguish stem cell enhancers (i. e. , non-enhancers at E11. 5) from those specific to differentiated lineages. We note that despite these interesting observations, most individual functional genomics features do not carry a great deal of information and the power of EnhancerFinder comes from the integration of different types of data. It is also possible that as a more complete experimental characterization of chromatin state and protein-DNA binding from E11. 5 tissues is obtained, data from less relevant contexts will not provide as much improvement as it did in this study. We focused on a single developmental stage with a large number of validated enhancers. To efficiently extend enhancer detection and validation to new contexts, it will be very important to select the most informative data to collect. Even though the ENCODE project has produced an impressive amount of data, it still has not extensively assayed most contexts of interest to researchers, in particular developmental biologists. The performance of classifiers trained on subsets of all our data and the weights we learned for feature sets and individual features provide some guidance for future experiments. Evolutionary conservation and DNA sequence patterns are broadly useful in the identification of enhancers, but our results suggest that adding functional genomics data is necessary to make more precise predictions about the contexts of activity. H3K4me1 and p300 are two of the most useful functional genomics data types overall (Figure S6), but many others are useful in particular contexts. However, the non-random sampling of functional genomics data and enhancers makes definitively determining the relative utility of different data types challenging. We saw a broad range in our ability to predict the tissue specificity of enhancers from existing data. Heart enhancers were dramatically easier to identify than other tissue-specific enhancers. Heart enhancers have significantly higher GC content than enhancers of other tissues, are less evolutionarily conserved, and are closer to the nearest TSS than other known enhancers at E11. 5, and we show that GC content alone is sufficient to accurately predict many heart enhancers (Figures S7Figures S7 and S8). However, functional genomics data alone were also able to accurately predict heart enhancers. The underlying biological explanation for these patterns may have to do with relative developmental age of different organs and tissues. At E11. 5, the heart is further along its developmental trajectory than the other tissues considered, and heart enhancers have completed their most conserved developmental stage, whereas forebrain enhancers are most strongly conserved at E11. 5 and E14. 5 [75]. At E11. 5, many of the less conserved, mammal-specific features of the heart are developing [91], [92], whereas other tissues are still developing under more general, less species-specific conserved regulatory programs at E11. 5 [93]. A recent study of enhancers in the adult mouse retina found that high local GC content was strongly correlated with enhancer activity [94]. Paired with our result, this suggests that GC content is a distinguishing feature of certain classes of enhancers. In spite of the strong overall performance of EnhancerFinder at predicting tissue-specific developmental enhancers, our approach has some limitations. First, we rely heavily upon the VISTA Enhancer Browser for training examples, because it is the largest collection of validated mammalian enhancers currently available. This resource provides an impressive catalog of validated human regulatory enhancers, but it is limited to a single developmental stage and experimental system. Without more data and analysis, it is difficult to evaluate how specific our predictions are to this context. Applying EnhancerFinder to known enhancers in model organisms, such as zebrafish and fly, would provide additional opportunities to evaluate our approach and findings, while potentially demonstrating differences in how enhancers function in these different species. Second, most of the enhancers present in VISTA are evolutionarily conserved. As a result, the VISTA enhancers cannot be viewed as an exhaustive catalog of the full range of enhancers. However, these regions have validated enhancer activity in vivo, and thus provide an appealing alternative to approaches that use single-mark proxies for enhancer activity (e. g. , considering all H3K27ac peaks as active enhancer regions). In addition to being conserved, these regions contain many signatures of enhancers in their sequence motifs and functional genomics composition that are useful for predicting enhancers. To emphasize these features and mitigate the impact of bias towards conserved regions, we removed evolutionary conservation as a feature from EnhancerFinder when we applied it to predict enhancers genome-wide. Our goal in doing so was to improve our ability to discern less conserved enhancers in these genome-wide predictions, and indeed, we predicted thousands of non-conserved enhancers (∼20% of all predictions). Third, though our predictions are based on a large collection of genome-wide chromatin state, protein-binding, and sequence information from many contexts, we are still limited by data availability. Even with the impressive efforts of ENCODE and related projects, producing data that are perfectly matched to all contexts of interest is time consuming and sometimes impossible, especially when studying humans. Thus, it will be important to develop a principled understanding of how different data can be generalized across tissues, developmental stages, and between species. In our analysis, many of the highest weighted features come from contexts close to the developmental stage of interest, and thus we anticipate that gathering more data from developmentally relevant cells and tissues will significantly improve our ability to annotate genomic regions involved in the regulation of embryonic development. However, data from other, seemingly unrelated, contexts may continue to prove useful. This study annotates regulatory elements in the human genome and provides tools for interpreting the effects of mutations in non-coding regions. Our case studies on regions around ZEB2, FOXC1, and FOXC2 illustrate how our predictions can facilitate the rapid identification of novel enhancers. In addition, the statistical enrichment for GWAS SNPs in our genome-wide enhancer predictions suggests that they may be a good resource for pinpointing causal mutations in potential disease loci. EnhancerFinder is a general framework for enhancer prediction and evaluation of different data sources that aim to annotate the regulatory functions of the human genome. It could easily be extended to include additional types of data, such as population-level variation at each locus, information about the three-dimensional state of the genome from Hi-C and 5C, and predictions of potential target genes for each enhancer. It could also be used to analyze additional aspects of the data we already consider, such as accounting for the relative genomic position of different features [66]. The EnhancerFinder two-step approach enables delineation of features common to all enhancers versus those that characterize enhancers of different types. For example, we find that predicting enhancers that are unique to a single tissue is more difficult than those that are active in multiple tissues (Figure S9), that certain features make prediction of heart enhancers particularly easy, and that different features are selected in classifiers for general enhancers and those for specific tissues. Together, these results suggest that there may be distinct classes of enhancers, even among those active in a given tissue at a single developmental stage. Further analysis of EnhancerFinder classifiers based on different types of data may help suggest biological mechanisms underlying the functional distinctions and genomic features of these different classes of enhancers.
Transgenic mice were generated by Cyagen Biosciences (http: //www. cyagen. com/). Their facility meets and often exceeds animal health and welfare guidelines. Animals were euthanized using techniques recommended by the American Veterinary Medical Association. All procedures were carried out in line with Gladstone Institutes and University of California guidelines. All zebrafish work was approved by the UCSF Institutional Animal Care and Use Committee (protocol number AN100466). All work presented in this paper is based on the February 2009 assembly of the human genome (GRCh37/hg19) downloaded from the UCSC Genome Browser (http: //genome. ucsc. edu/). Any data that was not in reference to this build was mapped over using the liftOver tool from the UCSC Kent tools (http: //hgdownload. cse. ucsc. edu/admin/jksrc. zip). In our framework, genomic regions are associated with a common set of descriptive features. We then apply machine-learning algorithms that use the features of known training examples to learn a function of the feature data that distinguishes the positives (enhancers) from the negatives (non-enhancers). This function can then be applied to the features associated with uncharacterized genomic regions to predict their enhancer status. A positive score for a genomic region indicates predicted membership in the positive class (enhancers) and a negative score indicates predicted membership in the negative class (non-enhancers). To evaluate the performance of trained classifiers, we performed 10-fold cross-validation on the training data and quantified our results with ROC AUC, precision-recall curves, and power estimates at fixed false positive rates. We computed p-values for the difference in performance between classification methods using McNemar' s test [99], [100]. To estimate false discovery rates, we trained EnhancerFinder classifiers at 1∶1,1∶10, and 1∶100 ratios of positive to negative enhancers and used the resulting 10-fold cross-validation results to calculate the proportion of false discoveries genome-wide at a 5% FPR if the true proportion of 1. 5 kb windows containing an enhancer was 50%, 10%, or 1%. We compared EnhancerFinder' s predictions to those of several previous enhancer prediction methods. We obtained the performance of CLARE on our Step 1 prediction task, by inputting our positive and negative data into the CLARE web server [73]. We downloaded the genomic segmentations and annotations produced by ChromHMM [64] and Segway [65]. We considered the ChromHMM predictions based on different ENCODE cell lines both individually and together. Any genomic region in our evaluation data set that overlapped an enhancer state was considered a predicted enhancer, and all others were considered predicted non-enhancers. For Segway, we also considered the “TF activity” state. We predicted tissue-specific developmental enhancers throughout the human genome by applying a trained MKL classifier (Step 1 of EnhancerFinder) without conservation (see Results) to sliding windows of 1500 bp, moving along the human genome in 500 bp steps. The feature profile for each window was computed as described above. To focus on high-confidence predictions, we filtered the enhancer scores for the windows at a 5% FPR, estimated from cross-validation using the genomic background, and combined the remaining overlapping windows to produce 84,301 high-confidence predicted enhancers. To predict tissue specificity, we applied trained brain, limb, and heart classifiers (Step 2 of EnhancerFinder) without conservation to all 299,039 windows with positive enhancer scores in Step 1. We then applied a 5% FPR cutoff for each tissue and concatenated the remaining overlapping windows into merged enhancer regions. Using this approach, we predicted 19,051 heart enhancers, 11,693 brain enhancers, and 7,400 limb enhancers. We characterized the expression patterns of the gene nearest to each predicted enhancer using the GNF Atlas 2 [101]. It contains expression data for genes in 79 different tissues, with expression measured using Affymetrix microarrays. For each of these 79 tissues, we used a paired t-test to determine if the nearest genes of predicted heart enhancers had significantly different mean values of expression than the nearest genes of brain enhancers. We did not include the limb enhancers in this analysis due to the lack of relevant expression data in the GNF Atlas 2. We examined genomic regions near predicted developmental enhancers for enrichment of Gene Ontology functional annotations, known phenotypes, and pathways using GREAT [102]. Results were computed using the hypergeometric test for genome-wide significance, with the default settings and the “basal plus extension” association rule (proximal 5 kb upstream, 1 kb downstream, plus distal up to 100 kb). We identified the sequence motifs present in each set of enhancers using the FIMO tool (Find Individual Motif Occurrences) from the MEME Suite of sequence motif analysis tools [103]. We considered known transcription factor binding motifs from the April 2011 release of the TRANSFAC database with a FIMO score threshold of 10e-5. We identified those occurrences that fell in predicted enhancers, and summarized motifs to identify the most prevalent TFs in each tissue-specific set of enhancers. We analyzed the overlap of predicted enhancers with GWAS SNPs, based on the NHGRI catalog of 9,687 GWAS SNPs downloaded from the UCSC Genome Browser in October 2012. Unadjusted permutation p-values were calculated by randomizing genomic locations of predicted enhancers (matching for length and chromosome, and avoiding assembly gaps) and overlapping these randomized regions with GWAS SNPs to assess significance of overlapping regions. Mouse enhancer assays were carried out in transient transgenic mouse embryos generated by pronuclear injections of enhancer assay constructs into FVB embryos (Cyagen Biosciences). Human and chimpanzee DNA sequences were inserted upstream of a minimal promoter Hsp68 and a LacZ reporter gene. The human sequence was amplified using primers 5′-TGTATGAAACCTGTTCACTCTCC-3′ and 5′-GCTTAAAACAACTACTAGAATCAGGC-3′ from the bacterial artificial chromosome (BAC) RP11-107E5 (from the BacPac resource at CHORI). The chimpanzee sequence was amplified using primers 5′-TGTATGAAACCTGTTCACTCTCC-3′ and 5′-GCTTAAAACAACTACTAGAATCAGGC-3′ from BAC CH251-677E03a (CHORI). The embryos were collected and stained for LacZ expression at E11. 5. Following the annotation policies of the VISTA Enhancer Browser, we required that consistent spatial expression patterns be present in three or more embryos with staining in order for the region to be considered an enhancer. Zebrafish enhancer assays were performed in transient transgenic zebrafish embryos. We tested candidate enhancer regions that ranged in length from 987 bp to 3,633 bp (see Table S6 for hg19 genomic coordinates), which we manually demarcated from within larger predicted enhancer regions based on signatures of likely enhancer function (including DnaseI hypersensitivity sites, transcription factor binding sites, histone modifications, and conservation). We performed PCR to obtain the candidate enhancer sequence using human genomic DNA (Roche). These were cloned into the E1b-GFP-Tol2 enhancer assay vector containing an E1b minimal promoter followed by GFP [104], and the construct was verified by sequencing. Each construct was injected with Tol2 mRNA into at least 100 single-cell fertilized zebrafish embryos. We annotated GFP expression at approximately 24 and 48 hours post fertilization (hpf), and considered an enhancer to be positive if we observed consistent expression in at least 15% of all fish alive at either 24 or 48 hpf [105], and suggestive of enhancer activity if we observed consistent expression in at least 10% of all fish alive at 24 or 48 hpf, after subtracting out percentages of tissue expression in fish injected with the empty enhancer vector. For each construct, at least 50 fish were analyzed for GFP expression at 48 hpf.
The human genome contains an immense amount of non-protein-coding DNA with unknown function. Some of this DNA regulates when, where, and at what levels genes are active during development. Enhancers, one type of regulatory element, are short stretches of DNA that can act as “switches” to turn a gene on or off at specific times in specific cells or tissues. Understanding where in the genome enhancers are located can provide insight into the genetic basis of development and disease. Enhancers are hard to identify, but clues about their locations are found in different types of data including DNA sequence, evolutionary history, and where proteins bind to DNA. Here, we introduce a new tool, called EnhancerFinder, which combines these data to predict the location and activity of enhancers during embryonic development. We trained EnhancerFinder on a large set of functionally validated human enhancers, and it proved to be very accurate. We used EnhancerFinder to predict tens of thousands of enhancers in the human genome and validated several of the predictions near three important developmental genes in mouse or zebrafish. EnhancerFinder' s predictions will be useful in understanding functional regions hidden in the vast amounts of human non-coding DNA.
developmental biology
genomics
functional genomics
organism development
genome analysis
gene expression
genetics
biology and life sciences
comparative genomics
molecular genetics
computational biology
Integrating Diverse Datasets Improves Developmental Enhancer Prediction