BACKGROUND
the hedgehog  signaling pathway is critical for embryonic organ development and adult tissue homeostasis across animal phyla .
in multiple tissue-specific settings, hh signaling directs specific cell fate choices, controls tissue patterning and governs cell proliferation.
in mammals, hh signaling originates with any of three ligands , indian hedgehog  or desert hedgehog ) .
release of the lipid modified hh ligand has been shown to be facilitated by scube in conjunction with the transmembrane protein dispatched .
once discharged, hh ligands interact with the patched  receptor protein and with hh-binding proteins boc, cdo and gas <dig> on target cells to relieve ptch-dependent inhibition of the smoothened  transmembrane protein  <cit> .
hh ligands can also be sequestered by the hedgehog-interacting protein , which dampens signaling  <cit> .
hh-mediated signal transduction culminates in the nucleus, with the binding of zinc-finger transcription factors  to target gene sequences  <cit> .
however, proteolytic processing determines whether the gli proteins act as repressors or activators.
gli <dig>  which is not processed, functions exclusively as a transcriptional activator and may act to amplify hh signals  <cit> .
gli <dig> and gli <dig> can be converted to a repressor form in the absence of hh ligand.
in the presence of the hh ligand, this processing is inhibited, allowing full-length gli proteins to traffic to the nucleus and activate gene expression .
processing of gli proteins requires passage through the cilia  <cit> ; the kinesin kif <dig> helps to properly construct the cilium and is enriched at the cilium tip, along with gli and sufu   <cit> .
the hh signaling pathway is regulated by both positive and negative feedback.
indeed, a number of hh pathway components, including boc, cdo, gas <dig>  gli <dig>  hhip, ptch <dig> and ptch <dig> are thought to be direct transcriptional targets of hh signaling in multiple tissue contexts .
thus, an important aspect of hh pathway self-regulation is integrated at the level of the enhancers that control response of the pathway target genes to local hh signaling levels.
however, despite the high functional conservation of this pathway, surprisingly little is known about the enhancer elements that control self-regulation in any organism.
one way to identify hh target enhancers is to perform chromatin immunoprecipitation .
genetically modified mouse models carrying inducible flag-tagged gli proteins have allowed analysis of gli binding sites in vivo in several different tissue contexts.
four in vivo gli binding studies, including three chip-chip analyses  <cit>  and one chip-seq study  <cit> , have been carried out using these models.
interestingly, examination of all four datasets for common gli binding sites that are annotated to hh pathway molecules reveals only three such sites  that are uniformly detectable.
several other established hh pathway genes, including boc, hhip, gli <dig>  and hipk <dig>  appear to exhibit different gli-bound genomic locations, depending on context, suggesting that each of these pathway components is regulated by multiple distinct genomic enhancers that have context-specific features.
using chip studies on diverse tissues, it may be possible to eventually identify all of the multiple enhancers that control each target gene in every context.
while a valuable goal, such analyses are currently expensive and time consuming and often technically challenging where the number of cells available for analysis is limiting, as in many developmental contexts.
importantly, computational methods can reveal sequence features that characterize enhancer activity.
we therefore asked whether analysis of all existing gli chip data could reveal common sequence features that might be used to empirically and globally predict functional enhancers de novo.
a publicly available machine learning approach, kmer-svm  <cit> , was used to predict novel hh enhancer regions.
this tool uses a support vector machine  to determine sequence features  that identify positive genomic regulatory regions  <cit> .
svms are classifier algorithms that define a boundary between members of two different groups.
kmer-svm calculates weights for sequence features that determine the effectiveness of that feature to distinguish between positive and negative regulatory regions.
once the features are determined, they can be used to identify novel enhancer regions not present in the original positive set.
the strength of this approach is that it relies exclusively on short regions of dna sequence  which are in the size range of transcription factor binding sites .
additionally, the organization of the k-mers within a sequence does not impact the score; this feature is consistent with the variable arrangement of tfbs in enhancers  <cit> .
using the kmer-svm tool  <cit> , analysis of the four existing gli binding datasets identified a set of k-mers that appeared to successfully predict potential gli-regulated enhancers.
application of this set of k-mers to the mouse genome pinpointed over  <dig>  potential enhancers.
several putative enhancers that were annotated to hh pathway components were then tested for their ability to drive hh-dependent activity in transfected cells.
the functional significance of the gli binding motifs  was also tested within each active enhancer by mutation.
of the nine predicted regulatory regions tested, seven  drove reporter expression in a gli-dependent fashion.
these findings substantially increase the number of functionally verified hh enhancers found in hh pathway molecules and validate the use of machine learning on chip data as a valuable tool to empirically predict likely hh-dependent regulatory regions.
RESULTS
analysis of gliflag datasets to identify likely in vivo gli transcription factor binding motifs
a previous in vitro analysis of gli transcription factor binding resulted in the identification of a set of likely binding sites for this factor  <cit> .
however, this spectrum of sites may not accurately represent gli binding site preferences in vivo.
to begin to examine this, we performed a meta-analysis of four existing gli-chip datasets.
all of these datasets utilize transgenic mice carrying flag tagged gli <dig>  or gli <dig>  in the rosa <dig> locus, activated by cre recombination, in four different tissue contexts: limb bud development   <cit> , cerebellum development , medulloblastoma  resulting from hh signaling overexpression  <cit> , and neural progenitor cells   <cit> .
an additional study of neural progenitors  <cit>  was excluded from analysis since it contained a low number of significant peaks and mirrors the same experimental conditions as the np dataset  <cit> .
for each of the datasets, the reported percentage of chip peak sequences with gli binding motifs  was as follows: ld 55 %, cd 26 %, mb 46 %, and np 91 %.
however, the definition of gbm was not the same across all datasets: one study allowed only two mismatches from the consensus  <cit>  and others generated a gli motif de novo based on the sequences of recovered peaks  <cit> .
to collate the spectrum of gbm observed in all four datasets, we applied a de novo motif enrichment analysis to each dataset individually  <cit> .
sequences that contained at least one site that matched the de novo motifs were removed from the dataset.
the remaining sequences were analyzed for residual motifs that resembled a gbm using dreme  <cit>  and tomtom  <cit>  .
this resulted in  <dig> putative gbm  , encompassing the range of gbm that are present in existing chip data.
this set therefore represents a collection of likely genomic gli binding sites, although some functional gli binding sites in vivo could be absent from this set and some false positive sites may be included.
each 12-mer was classified as high confidence , medium confidence , or low confidence  if it was found within sequences from all four datasets, two to three datasets, or one dataset, respectively.
the sequence logos  <cit>  for each classification, provided in fig.
1a, show a nearly absolute representation of ccxc in positions 4– <dig> for all sites.
indeed, concordant  nucleotides at the 5th and 7th position were previously found to be required for gli binding  <cit> .
interestingly, for high confidence sites, there is no variation at  <dig> of the  <dig> positions, including the 5th and 7th positions .fig.
<dig> definition of gli binding motifs  and characterization of gliflag datastets.
a sequence logos  of 12-mer gbm.
high confidence , medium confidence , and low confidence  gbm are found in all four , three or two  or one  datasets.
b relative frequency of peaks containing one or more gbm in the gliflag chip-chip ld , cd , mb , and gliflag chip-seq np  datasets.
a high proportion of sequences contain only one gbm.
c overlap of sequences identified by all four gliflag datasets; only  <dig> individual peaks are found in all contexts
using the recommended matrix similarity score cutoff of 81 % overall matrix similarity to the optimal consensus gli site as defined by the in vitro dna binding assay  <cit>  results in  <dig> , <dig> putative gli tfbs across the mouse genome.
this is substantially more than the  <dig>  found using the new gli library defined by the in vivo chip studies.
however, several of the sites predicted by the in vitro binding studies do not contain the concordant  nucleotides at the 5th and 7th position.
thus, the newly generated gli library  may more accurately represent functional gbm.
according to this new library of  <dig> gbm, 41 % of ld, 27 % of cd, 32 % of mb, and 80 % of np peaks contain putative gli binding sites.
for those sequences that contain a gbm, the vast majority contain only a single site  .
the overlap of genomic binding regions among datasets is shown in fig.
1c; only  <dig> genomic coordinates are shared among all datasets .
since pathway components must respond to the hh signal in all tissues, it might be expected that this common response would be integrated by a single enhancer.
however, only three of the  <dig> shared regions are annotated to known hh pathway components .
indeed, for boc and hhip, distinct genomic gli binding regions are found in different datasets.
this result suggests that some of these genes may have multiple enhancers that work to transduce the hh signal in different tissue contexts.
given this apparent complexity in regulatory regions, we next asked whether the existing datasets of chip peaks might contain additional sequence information that could be used to predict the location of other hh-responsive enhancers in the mouse genome.
a machine learning approach was employed to test this question.
assessment of kmer-svm performance and prediction
kmer-svm assessment of classification using gliflag datasets
for each gliflag dataset, only sequences with at least one gbm  were used.
this was done since a high proportion of the chip-chip datasets did not contain a putative gli binding site as defined in the original papers  or by our assessment .
each individual dataset was submitted to kmer-svm and the ability of each classifier to correctly label a candidate sequence as positive was assessed.
background sequences were randomly selected from the genome, but matched for gc content with the positive set.
kmer-svm randomly divides the data as follows: 80 % of the sequences are used as a training set and 20 % are used as a testing set.
the ability of the classifier built with the training set to accurately identify the members of the remaining 20 % testing set is then assessed.
this is repeated five times, each with a different random division of the data.
receiver operating characteristic  curves and precision recall curves  are used to assess the success of the classifier to correctly label regions in the testing set as positive .
roc curves display the cumulative distribution of the true positive rate compared to the false positive rate.
this characteristic assesses how well the classifier is able to label the positive sequences from the test set.
the area under the curve was  <dig>  for ldwgbm ,  <dig>  for cdwgbm ,  <dig>  for mbwgbm  and  <dig>  for npwgbm .
thus, the classifier performs best in ld and np datasets.fig.
<dig> assessment of classification capability of kmer-svm trained gliflag datasets containing sequences with at least one gbm.
for all curves, each dataset is randomly split into 80 % for training and 20 % for prediction and the prediction is repeated five times .
plots assess the likelihood that the specified classifier can successfully predict sequences that have at least one gbm as positive or negative.
a-d roc plots depicting true positive rates  and false positive rates .
area under the curve  scores as calculated by kmer-svm are:  <dig>  for ldwgbm ,  <dig>  for cdwgbm ,  <dig>  for mbwgbm  and  <dig>  for npwgbm  datasets.
e-h precision recall curves depicting the positive predictive value , calculated as true positive / , versus the tpr.
auc of  <dig>  for ldwgbm  and  <dig>  for npwgbm  indicate reasonable confidence in the classification while auc of  <dig>  for cdwgbm  and  <dig>  for mbwgbm  indicate a low probability that the region is correctly labeled when the sequence is classified as positive
prc displays the predictive value against the true positive rate and represents the accuracy of the labeling.
the prc plots indicated high values for ldwgbm  and npwgbm  but low values for cdwgbm   and mbwgbm  .
the roc and prc plots for ldwgbm and npwgbm datasets suggested that the classifier sequence features used were able to distinguish between positive and negative groups with a low level of false labeling.
predictions
the classifiers for ldwgbm and npwgbm were then individually run on 600 bp of sequence centered on every gbm in the mouse genome .
use of both the ldwgbm and npwgbm datasets for prediction incorporated data from the gli1flag  and gli3flag  transcription factors in two diverse contexts .
the length of 600 bp was selected based on motif enrichment analysis of the ld and np datasets using meme-chip  <cit>  and centrimo  <cit> .
this analysis showed that, within the chip-chip ld dataset, enrichment for the location of gli motifs  has a broad profile that spans 200 bp to either side of the midpoint .
the gli motif has a narrower profile in the np data, a feature that is expected for chip-seq .
the profile for the sox motif , an established tissue specific gli cofactor  <cit> , shows an enrichment peak that is centered around 200 bp on either side of the midpoint  and suggests that cofactors for hh may reside outside of the immediate vicinity of a gli binding site.
we therefore used 600 bp to capture both common hh features as well as potential context specific sequence.
for the ldwgbm classifier, scores ranged from − <dig>  to  <dig>  with  <dig>  % of the  <dig>  analyzed genomic regions scoring as positive .
the npwgbm results ranged from − <dig>  to  <dig>  with  <dig>  % positive .
the categorization of a sequence is dependent on the sign of the score and the weight of the value is less important than the ranking.
overall, the correlation between scores for individual genomic regions calculated by the ldwgbm and npwgbm classifiers is poor  .
however, if only sequences with positive scores are considered, the correlation improves .
if scores are restricted to values indicating only the high confidence scores , the values are very well correlated .
in total,  <dig> genomic regions were predicted as hh enhancers by both classifiers.
of those,  <dig> regions  overlapped at least one peak in the four gliflag datasets.
among high confidence scores   <dig> genomic regions were shared between the two datasets and  <dig> of these  overlapped with peaks from at least one of the four gliflag datasets.
all of the scored regions are listed in additional file 4: table s <dig> fig.
<dig> assessment of genomic kmer-svm predictions using classifiers trained on ldwgbm and npwgbm datasets.
a all genomic sequences matching the restricted  <dig> gbm 12-mers  were identified and the 600 bp surrounding each gbm were assessed and scored using the kmer-svm classifier that was trained on each of the two datasets.
b correlation plot depicting the relationship between ldwgbm and npwgbm scores; scores > <dig> are highly correlated in the two datasets.
c gli motif generated from overlapping high weighted k-mers shared between ldwgbm and npwgbm classifiers.
d high weighted k-mers  represented in either ldwgbm  or npwgbm
evaluation of predictions
to assess whether kmer-svm predictions were likely to represent hh enhancers, we examined whether the predicted genomic regions overlapped publically available open chromatin and enhancer histone marks in tissues that were likely to be transducing hh signals.
we first examined the dnasei hypersensitive profile collected from mouse mesoderm at e <dig>  at genomic regions that were predicted with high confidence as positive  or negative  .
a significantly higher proportion of overlap was found with the predicted positive regions than with predicted negative regions  .
we also examined publically available h3k4me <dig>  h3k4me <dig>  and h3k27ac encode data collected from heart and liver at e <dig> .
monomethylated h3k <dig>  and histone h <dig> acetyl lys <dig> h3k27ac  <cit>  were used as enhancer markers while trimethylated h3k <dig>  was expected to be depleted in enhancer regions  <cit> .
although hh signaling is active during early development of both tissues, available in situ analysis for gli <dig>   <cit>  shows gli <dig> expression in liver but not heart at e <dig>  .
thus, we expected signals for both h3k4me <dig>  and h3k27ac  to be enriched in the predicted positive regions in liver but not heart at this time point.
indeed this was the case: h3k4me <dig>  and h3k27ac , with no significant difference in h3k4me <dig> when predicted positive regions were compared to predicted negative regions.
as expected, the heart data did not show enrichment for h3k4me <dig> or h3k27ac.
together, the results  are consistent with the conclusion that the kmer-svm classification correctly identifies hh enhancer regions.
next, we evaluated the sequence features, or k-mers, that kmer-svm identified as primary components of hh enhancer regions.
the weights of k-mers are calculated during the svm training and reflect the contribution of the k-mer to categorization of a sequence.
weights can be positive or negative and the sum of the weights of iterative k-mers across a sequence comprise the overall score of that sequence.
not surprisingly, alignment of k-mers with high scoring weights shared between both datasets returned a motif that strongly resembles the gbm .
unique high weighted k-mers that occurred in each individual dataset represented potential context specific features.
an e-box motif was identified for the ldwgbm dataset while a sox motif was returned for npwgbm .
negative weights that occurred in both datasets include ac and acc repeats as well as other c rich sequences.
functional verification of gli-dependent enhancer activity
predicted genomic regions were annotated to the two nearest genes using great  <cit> .
because our goal was to identify enhancers for hh pathway components, we selected a subset of predictions that were positive in both the ldwgbm and npwgbm datasets and that were annotated to members of the go: <dig> smoothened signaling pathway gene set.
because hh pathway components are required for active hh signaling, we reasoned that enhancers annotated to these genes would be more likely to function in any tissue that transduces hh signal.
therefore, high scoring regions annotated to different members of the go: <dig> gene set that were readily cloned were functionally tested for enhancer activity.
two previously known hh enhancers for ptch <dig> and ptch <dig> appeared on this list  <cit> .
interestingly, an established gli <dig> regulatory region was not predicted  <cit> .
the test set consisted of genomic regions annotated to boc, gli <dig>  hhip, hipk <dig>  ptch <dig>  scube <dig>  shh, and tgfbr <dig>  an additional region, annotated to dpp <dig>  was also tested .table  <dig> assessment of predicted hh enhancer regions
ptch2
boc
dpp6
gli3
hhip
hipk2
ptch1
scube1
shh
tgfbr2
seven of the nine regions predicted to be gli-driven enhancers were indeed determined to be hh responsive and gli binding site dependent in a cell culture assay.
overlap of the predicted regions with peaks from the gliflag chip datasets  is indicated by the plus sign.
boc, hipk <dig>  scube <dig> and tgfbr <dig> were predicted by kmer-svm and found to be positive, even though those regions do not overlap any of the peaks identified in the ld or np datasets used to generate the kmer-svm classifier
the nine genomic regions were screened for hh responsiveness in a cell culture assay that has been previously used to validate hh enhancers  <cit> .
in this assay, c2c <dig> cells are transiently transfected with constructs containing the putative hh regulatory region cloned upstream of a minimal promoter that drives luciferase expression .
to induce a hh response, enhancer reporter constructs were co-transfected with a plasmid that drives constitutive expression of gli <dig> in mammalian cells  <cit> .
for those candidate enhancers that demonstrated apparent hh activation, gli-dependent activity was confirmed by retesting after mutagenesis of the gli binding sites .
either complete loss of enhancer activity or attenuation of response in gko sequence was considered gli-dependent.
the established hh enhancer region for ptch <dig> was used as a positive control  <cit> .
for the  <dig> regulatory regions annotated to hh pathway component genes,  <dig> exhibited hh activation that was directly dependent on a gli binding site .
the ptch <dig> positive control region showed a complete loss of hh response after mutation of the gli site  as did the regulatory regions annotated to hhip, hipk <dig>  ptch <dig>  and scube <dig>  regions annotated to boc, dpp <dig>  and tgfbr <dig> showed a significant decrease in hh activation upon gli tfbs mutation, but not a complete loss of response.
this suggests that additional regulatory inputs influence the activity of these enhancers.
neither the gli <dig> nor shh region exhibited hh dependent enhancer activity.
however, we cannot rule out the possibility that these regions might be positive if examined in a different cellular context  <cit> .fig.
<dig> functional verification of gli-dependent enhancer activity.
putative regulatory regions were cloned upstream of a minimal promoter co-transfected into c2c <dig> cells, along with a gli <dig> expression vector.
relative activity is plotted .
the ptch <dig> region is used as a positive control.
novel regions annotated to the boc, dpp <dig>  hhip, hipk <dig>  ptch <dig>  scube <dig>  and tgfbr <dig> loci exhibit upregulation in response to gli <dig> co-transfection .
dependence on gli was assessed by mutating all putative gli tfbs  within the sequence and retesting in the assay
two previous in silico methods have been described for the identification of hh-regulated enhancers in vertebrates: enhancer element locator  and module cluster analysis .
eel analyzes the estimated energy of a single transcription factor binding event, as well as the possible interaction between adjacent, highly conserved transcription factor binding sites, to detect potential enhancers  <cit> .
in contrast, mca utilizes poisson modeling to determine the relative enrichment of binding sites in highly conserved, non-coding sequence and, thereby, identify putative enhancers  <cit> .
while both of these methods have had some success  in detecting putative hh-regulated enhancers, a disadvantage to these approaches is that the search is limited to regions of high sequence conservation and to regions close to promoters.
in contrast, kmer-svm approach used here employs a genome-wide empirical analysis to locate regions that contain sequence features predictive of hh enhancer function.
though these predictions miss one of the three known hh pathway enhancers  <cit>  indicating that the algorithm does not capture all hh-driven enhancers, the high success rate  of the kmer-svm predictions far exceed the previous prediction rates for eel or mca.
it is important to note, however, that the predictions tested here are all for pathway genes, which may have a unique signature.
it would be necessary to test additional enhancers in tissue-specific assay systems  to determine the overall success of this method in identification of tissue-specific enhancers.
to examine the impact of k-mers that contribute to predictions of the positive regions, weights were plotted across each of the sequences.
the ptch <dig> sequence , a known enhancer region  <cit> , contained matching profiles for 8-mers predicted from ldwgbm  and npwgbm  that are gc rich and similar to the gbm k-mers however, mutation of the single gli tfbs  ablates the hh response, indicating that the presence of this gbm is required to transduce hh signaling.
mutation of the gbms annotated by the green boxes for hhip, hipk <dig>  ptch <dig>  and scube <dig>  is also sufficient to abrogate hh signal transduction.
for boc, dpp <dig>  and tgfbr <dig>  which show enhancer activity that remains after ablation of the gbm , there were no sequence characteristics that were indicative of a shared feature responsible for this remaining response to induction by gli <dig>  the boc profile was the only one that contained a high weighted k-mer  that was unique to the ldwgbm  profile.
this k-mer was similar to a krox motif   <cit>  and may be enriched in the ldwgbm dataset as a context specific transcription factor, since it has roles in limb development  <cit> .
in general, most of the tested regions contain distributed high weighted 8-mers in addition to the central gbm and had profiles that contained consistent peaks in both ldwgm and npwgbm datasets .fig.
<dig> k-mer weights plotted across sequences that show enhancer activity.
diagrams were generated in ucsc genome browser and show coordinate information for regions annotated to ptch <dig> , boc , dpp <dig> , gli <dig> , hhip , hipk <dig> , ptch <dig> , scube <dig> , shh  and tgfbr <dig> .
green boxes represent gbm.
weights for ldwgbm and npwgbm are represented by the red and yellow lines, respectively.
refseq gene annotations are represented in blue.
a putative krox- <dig> tfbs  that has a high weight in the ldwgbm classifier but not the npwgbm classifier occurs in the sequence annotated to boc.
note that most sequences show weighted k-mers located several hundred bp from the central gbm, suggesting that sequence motifs that predict hh enhancer activity may be distributed
CONCLUSIONS
we have utilized the kmer-svm machine learning approach to examine four existing gli chip databases and to generate classifiers that can empirically predict functional hh enhancers from genomic dna.
the analysis was facilitated by a new gbs library generated from a meta-analysis of genomic regions identified in in vivo binding studies .
when compared to the previous library derived from in vitro binding studies  <cit> , this new gbs library resulted in the identification of a subset  of potential gli transcription factor binding sites across the mouse genome.
of nine predicted hh target regulatory regions tested, seven were functionally verified as gli-dependent.
all of the tested regions were annotated to genes believed to be components of the hh pathway and important determinants of the hh response.
given the high success rate of hh enhancer prediction in this small scale screen , it is quite possible that a large number of the other  <dig>  predicted regions  may harbor hh enhancer activity.

 the hedgehog signaling pathway is critical for embryonic organ development and adult tissue homeostasis across animal phyla. in multiple tissue-specific settings, hh signaling directs specific cell fate choices, controls tissue patterning and governs cell proliferation. in mammals, hh signaling originates with any of three ligands, indian hedgehog or desert hedgehog ). release of the lipid modified hh ligand has been shown to be facilitated by scube in conjunction with the transmembrane protein dispatched. once discharged, hh ligands interact with the patched receptor protein and with hh-binding proteins boc, cdo and gas on target cells to relieve ptch-dependent inhibition of the smoothened transmembrane protein. hh ligands can also be sequestered by the hedgehog-interacting protein, which dampens signaling. hh-mediated signal transduction culminates in the nucleus, with the binding of zinc-finger transcription factors to target gene sequences. however, proteolytic processing determines whether the gli proteins act as repressors or activators. gli which is not processed, functions exclusively as a transcriptional activator and may act to amplify hh signals. gli and gli can be converted to a repressor form in the absence of hh ligand. in the presence of the hh ligand, this processing is inhibited, allowing full-length gli proteins to traffic to the nucleus and activate gene expression. processing of gli proteins requires passage through the cilia; the kinesin kif helps to properly construct the cilium and is enriched at the cilium tip, along with gli and sufu. the hh signaling pathway is regulated by both positive and negative feedback. indeed, a number of hh pathway components, including boc, cdo, gas gli hhip, ptch and ptch are thought to be direct transcriptional targets of hh signaling in multiple tissue contexts. thus, an important aspect of hh pathway self-regulation is integrated at the level of the enhancers that control response of the pathway target genes to local hh signaling levels. however, despite the high functional conservation of this pathway, surprisingly little is known about the enhancer elements that control self-regulation in any organism. one way to identify hh target enhancers is to perform chromatin immunoprecipitation. genetically modified mouse models carrying inducible flag-tagged gli proteins have allowed analysis of gli binding sites in vivo in several different tissue contexts. four in vivo gli binding studies, including three chip-chip analyses and one chip-seq study, have been carried out using these models. interestingly, examination of all four datasets for common gli binding sites that are annotated to hh pathway molecules reveals only three such sites that are uniformly detectable. several other established hh pathway genes, including boc, hhip, gli and hipk appear to exhibit different gli-bound genomic locations, depending on context, suggesting that each of these pathway components is regulated by multiple distinct genomic enhancers that have context-specific features. using chip studies on diverse tissues, it may be possible to eventually identify all of the multiple enhancers that control each target gene in every context. while a valuable goal, such analyses are currently expensive and time consuming and often technically challenging where the number of cells available for analysis is limiting, as in many developmental contexts. importantly, computational methods can reveal sequence features that characterize enhancer activity. we therefore asked whether analysis of all existing gli chip data could reveal common sequence features that might be used to empirically and globally predict functional enhancers de novo. a publicly available machine learning approach, kmer-svm, was used to predict novel hh enhancer regions. this tool uses a support vector machine to determine sequence features that identify positive genomic regulatory regions. svms are classifier algorithms that define a boundary between members of two different groups. kmer-svm calculates weights for sequence features that determine the effectiveness of that feature to distinguish between positive and negative regulatory regions. once the features are determined, they can be used to identify novel enhancer regions not present in the original positive set. the strength of this approach is that it relies exclusively on short regions of dna sequence which are in the size range of transcription factor binding sites. additionally, the organization of the k-mers within a sequence does not impact the score; this feature is consistent with the variable arrangement of tfbs in enhancers. using the kmer-svm tool, analysis of the four existing gli binding datasets identified a set of k-mers that appeared to successfully predict potential gli-regulated enhancers. application of this set of k-mers to the mouse genome pinpointed over potential enhancers. several putative enhancers that were annotated to hh pathway components were then tested for their ability to drive hh-dependent activity in transfected cells. the functional significance of the gli binding motifs was also tested within each active enhancer by mutation. of the nine predicted regulatory regions tested, seven drove reporter expression in a gli-dependent fashion. these findings substantially increase the number of functionally verified hh enhancers found in hh pathway molecules and validate the use of machine learning on chip data as a valuable tool to empirically predict likely hh-dependent regulatory regions.  analysis of gliflag datasets to identify likely in vivo gli transcription factor binding motifs a previous in vitro analysis of gli transcription factor binding resulted in the identification of a set of likely binding sites for this factor. however, this spectrum of sites may not accurately represent gli binding site preferences in vivo. to begin to examine this, we performed a meta-analysis of four existing gli-chip datasets. all of these datasets utilize transgenic mice carrying flag tagged gli or gli in the rosa locus, activated by cre recombination, in four different tissue contexts: limb bud development, cerebellum development, medulloblastoma resulting from hh signaling overexpression, and neural progenitor cells. an additional study of neural progenitors was excluded from analysis since it contained a low number of significant peaks and mirrors the same experimental conditions as the np dataset. for each of the datasets, the reported percentage of chip peak sequences with gli binding motifs was as follows: ld 55 %, cd 26 %, mb 46 %, and np 91 %. however, the definition of gbm was not the same across all datasets: one study allowed only two mismatches from the consensus and others generated a gli motif de novo based on the sequences of recovered peaks. to collate the spectrum of gbm observed in all four datasets, we applied a de novo motif enrichment analysis to each dataset individually. sequences that contained at least one site that matched the de novo motifs were removed from the dataset. the remaining sequences were analyzed for residual motifs that resembled a gbm using dreme and tomtom. this resulted in putative gbm, encompassing the range of gbm that are present in existing chip data. this set therefore represents a collection of likely genomic gli binding sites, although some functional gli binding sites in vivo could be absent from this set and some false positive sites may be included. each 12-mer was classified as high confidence, medium confidence, or low confidence if it was found within sequences from all four datasets, two to three datasets, or one dataset, respectively. the sequence logos for each classification, provided in fig. 1a, show a nearly absolute representation of ccxc in positions 4– for all sites. indeed, concordant nucleotides at the 5th and 7th position were previously found to be required for gli binding. interestingly, for high confidence sites, there is no variation at of the positions, including the 5th and 7th positions.fig. definition of gli binding motifs and characterization of gliflag datastets. a sequence logos of 12-mer gbm. high confidence, medium confidence, and low confidence gbm are found in all four, three or two or one datasets. b relative frequency of peaks containing one or more gbm in the gliflag chip-chip ld, cd, mb, and gliflag chip-seq np datasets. a high proportion of sequences contain only one gbm. c overlap of sequences identified by all four gliflag datasets; only individual peaks are found in all contexts using the recommended matrix similarity score cutoff of 81 % overall matrix similarity to the optimal consensus gli site as defined by the in vitro dna binding assay results in, putative gli tfbs across the mouse genome. this is substantially more than the found using the new gli library defined by the in vivo chip studies. however, several of the sites predicted by the in vitro binding studies do not contain the concordant nucleotides at the 5th and 7th position. thus, the newly generated gli library may more accurately represent functional gbm. according to this new library of gbm, 41 % of ld, 27 % of cd, 32 % of mb, and 80 % of np peaks contain putative gli binding sites. for those sequences that contain a gbm, the vast majority contain only a single site. the overlap of genomic binding regions among datasets is shown in fig. 1c; only genomic coordinates are shared among all datasets. since pathway components must respond to the hh signal in all tissues, it might be expected that this common response would be integrated by a single enhancer. however, only three of the shared regions are annotated to known hh pathway components. indeed, for boc and hhip, distinct genomic gli binding regions are found in different datasets. this result suggests that some of these genes may have multiple enhancers that work to transduce the hh signal in different tissue contexts. given this apparent complexity in regulatory regions, we next asked whether the existing datasets of chip peaks might contain additional sequence information that could be used to predict the location of other hh-responsive enhancers in the mouse genome. a machine learning approach was employed to test this question. assessment of kmer-svm performance and prediction kmer-svm assessment of classification using gliflag datasets for each gliflag dataset, only sequences with at least one gbm were used. this was done since a high proportion of the chip-chip datasets did not contain a putative gli binding site as defined in the original papers or by our assessment. each individual dataset was submitted to kmer-svm and the ability of each classifier to correctly label a candidate sequence as positive was assessed. background sequences were randomly selected from the genome, but matched for gc content with the positive set. kmer-svm randomly divides the data as follows: 80 % of the sequences are used as a training set and 20 % are used as a testing set. the ability of the classifier built with the training set to accurately identify the members of the remaining 20 % testing set is then assessed. this is repeated five times, each with a different random division of the data. receiver operating characteristic curves and precision recall curves are used to assess the success of the classifier to correctly label regions in the testing set as positive. roc curves display the cumulative distribution of the true positive rate compared to the false positive rate. this characteristic assesses how well the classifier is able to label the positive sequences from the test set. the area under the curve was for ldwgbm, for cdwgbm, for mbwgbm and for npwgbm. thus, the classifier performs best in ld and np datasets.fig. assessment of classification capability of kmer-svm trained gliflag datasets containing sequences with at least one gbm. for all curves, each dataset is randomly split into 80 % for training and 20 % for prediction and the prediction is repeated five times. plots assess the likelihood that the specified classifier can successfully predict sequences that have at least one gbm as positive or negative. a-d roc plots depicting true positive rates and false positive rates. area under the curve scores as calculated by kmer-svm are: for ldwgbm, for cdwgbm, for mbwgbm and for npwgbm datasets. e-h precision recall curves depicting the positive predictive value, calculated as true positive /, versus the tpr. auc of for ldwgbm and for npwgbm indicate reasonable confidence in the classification while auc of for cdwgbm and for mbwgbm indicate a low probability that the region is correctly labeled when the sequence is classified as positive prc displays the predictive value against the true positive rate and represents the accuracy of the labeling. the prc plots indicated high values for ldwgbm and npwgbm but low values for cdwgbm and mbwgbm. the roc and prc plots for ldwgbm and npwgbm datasets suggested that the classifier sequence features used were able to distinguish between positive and negative groups with a low level of false labeling. predictions the classifiers for ldwgbm and npwgbm were then individually run on 600 bp of sequence centered on every gbm in the mouse genome. use of both the ldwgbm and npwgbm datasets for prediction incorporated data from the gli1flag and gli3flag transcription factors in two diverse contexts. the length of 600 bp was selected based on motif enrichment analysis of the ld and np datasets using meme-chip and centrimo. this analysis showed that, within the chip-chip ld dataset, enrichment for the location of gli motifs has a broad profile that spans 200 bp to either side of the midpoint. the gli motif has a narrower profile in the np data, a feature that is expected for chip-seq. the profile for the sox motif, an established tissue specific gli cofactor, shows an enrichment peak that is centered around 200 bp on either side of the midpoint and suggests that cofactors for hh may reside outside of the immediate vicinity of a gli binding site. we therefore used 600 bp to capture both common hh features as well as potential context specific sequence. for the ldwgbm classifier, scores ranged from − to with  % of the analyzed genomic regions scoring as positive. the npwgbm results ranged from − to with  % positive. the categorization of a sequence is dependent on the sign of the score and the weight of the value is less important than the ranking. overall, the correlation between scores for individual genomic regions calculated by the ldwgbm and npwgbm classifiers is poor. however, if only sequences with positive scores are considered, the correlation improves. if scores are restricted to values indicating only the high confidence scores, the values are very well correlated. in total, genomic regions were predicted as hh enhancers by both classifiers. of those, regions overlapped at least one peak in the four gliflag datasets. among high confidence scores genomic regions were shared between the two datasets and of these overlapped with peaks from at least one of the four gliflag datasets. all of the scored regions are listed in additional file 4: table s fig. assessment of genomic kmer-svm predictions using classifiers trained on ldwgbm and npwgbm datasets. a all genomic sequences matching the restricted gbm 12-mers were identified and the 600 bp surrounding each gbm were assessed and scored using the kmer-svm classifier that was trained on each of the two datasets. b correlation plot depicting the relationship between ldwgbm and npwgbm scores; scores > are highly correlated in the two datasets. c gli motif generated from overlapping high weighted k-mers shared between ldwgbm and npwgbm classifiers. d high weighted k-mers represented in either ldwgbm or npwgbm evaluation of predictions to assess whether kmer-svm predictions were likely to represent hh enhancers, we examined whether the predicted genomic regions overlapped publically available open chromatin and enhancer histone marks in tissues that were likely to be transducing hh signals. we first examined the dnasei hypersensitive profile collected from mouse mesoderm at e at genomic regions that were predicted with high confidence as positive or negative. a significantly higher proportion of overlap was found with the predicted positive regions than with predicted negative regions. we also examined publically available h3k4me h3k4me and h3k27ac encode data collected from heart and liver at e. monomethylated h3k and histone h acetyl lys h3k27ac were used as enhancer markers while trimethylated h3k was expected to be depleted in enhancer regions. although hh signaling is active during early development of both tissues, available in situ analysis for gli shows gli expression in liver but not heart at e. thus, we expected signals for both h3k4me and h3k27ac to be enriched in the predicted positive regions in liver but not heart at this time point. indeed this was the case: h3k4me and h3k27ac, with no significant difference in h3k4me when predicted positive regions were compared to predicted negative regions. as expected, the heart data did not show enrichment for h3k4me or h3k27ac. together, the results are consistent with the conclusion that the kmer-svm classification correctly identifies hh enhancer regions. next, we evaluated the sequence features, or k-mers, that kmer-svm identified as primary components of hh enhancer regions. the weights of k-mers are calculated during the svm training and reflect the contribution of the k-mer to categorization of a sequence. weights can be positive or negative and the sum of the weights of iterative k-mers across a sequence comprise the overall score of that sequence. not surprisingly, alignment of k-mers with high scoring weights shared between both datasets returned a motif that strongly resembles the gbm. unique high weighted k-mers that occurred in each individual dataset represented potential context specific features. an e-box motif was identified for the ldwgbm dataset while a sox motif was returned for npwgbm. negative weights that occurred in both datasets include ac and acc repeats as well as other c rich sequences. functional verification of gli-dependent enhancer activity predicted genomic regions were annotated to the two nearest genes using great. because our goal was to identify enhancers for hh pathway components, we selected a subset of predictions that were positive in both the ldwgbm and npwgbm datasets and that were annotated to members of the go: smoothened signaling pathway gene set. because hh pathway components are required for active hh signaling, we reasoned that enhancers annotated to these genes would be more likely to function in any tissue that transduces hh signal. therefore, high scoring regions annotated to different members of the go: gene set that were readily cloned were functionally tested for enhancer activity. two previously known hh enhancers for ptch and ptch appeared on this list. interestingly, an established gli regulatory region was not predicted. the test set consisted of genomic regions annotated to boc, gli hhip, hipk ptch scube shh, and tgfbr an additional region, annotated to dpp was also tested.table assessment of predicted hh enhancer regions ptch2 boc dpp6 gli3 hhip hipk2 ptch1 scube1 shh tgfbr2 seven of the nine regions predicted to be gli-driven enhancers were indeed determined to be hh responsive and gli binding site dependent in a cell culture assay. overlap of the predicted regions with peaks from the gliflag chip datasets is indicated by the plus sign. boc, hipk scube and tgfbr were predicted by kmer-svm and found to be positive, even though those regions do not overlap any of the peaks identified in the ld or np datasets used to generate the kmer-svm classifier the nine genomic regions were screened for hh responsiveness in a cell culture assay that has been previously used to validate hh enhancers. in this assay, c2c cells are transiently transfected with constructs containing the putative hh regulatory region cloned upstream of a minimal promoter that drives luciferase expression. to induce a hh response, enhancer reporter constructs were co-transfected with a plasmid that drives constitutive expression of gli in mammalian cells. for those candidate enhancers that demonstrated apparent hh activation, gli-dependent activity was confirmed by retesting after mutagenesis of the gli binding sites. either complete loss of enhancer activity or attenuation of response in gko sequence was considered gli-dependent. the established hh enhancer region for ptch was used as a positive control. for the regulatory regions annotated to hh pathway component genes, exhibited hh activation that was directly dependent on a gli binding site. the ptch positive control region showed a complete loss of hh response after mutation of the gli site as did the regulatory regions annotated to hhip, hipk ptch and scube regions annotated to boc, dpp and tgfbr showed a significant decrease in hh activation upon gli tfbs mutation, but not a complete loss of response. this suggests that additional regulatory inputs influence the activity of these enhancers. neither the gli nor shh region exhibited hh dependent enhancer activity. however, we cannot rule out the possibility that these regions might be positive if examined in a different cellular context.fig. functional verification of gli-dependent enhancer activity. putative regulatory regions were cloned upstream of a minimal promoter co-transfected into c2c cells, along with a gli expression vector. relative activity is plotted. the ptch region is used as a positive control. novel regions annotated to the boc, dpp hhip, hipk ptch scube and tgfbr loci exhibit upregulation in response to gli co-transfection. dependence on gli was assessed by mutating all putative gli tfbs within the sequence and retesting in the assay two previous in silico methods have been described for the identification of hh-regulated enhancers in vertebrates: enhancer element locator and module cluster analysis. eel analyzes the estimated energy of a single transcription factor binding event, as well as the possible interaction between adjacent, highly conserved transcription factor binding sites, to detect potential enhancers. in contrast, mca utilizes poisson modeling to determine the relative enrichment of binding sites in highly conserved, non-coding sequence and, thereby, identify putative enhancers. while both of these methods have had some success in detecting putative hh-regulated enhancers, a disadvantage to these approaches is that the search is limited to regions of high sequence conservation and to regions close to promoters. in contrast, kmer-svm approach used here employs a genome-wide empirical analysis to locate regions that contain sequence features predictive of hh enhancer function. though these predictions miss one of the three known hh pathway enhancers indicating that the algorithm does not capture all hh-driven enhancers, the high success rate of the kmer-svm predictions far exceed the previous prediction rates for eel or mca. it is important to note, however, that the predictions tested here are all for pathway genes, which may have a unique signature. it would be necessary to test additional enhancers in tissue-specific assay systems to determine the overall success of this method in identification of tissue-specific enhancers. to examine the impact of k-mers that contribute to predictions of the positive regions, weights were plotted across each of the sequences. the ptch sequence, a known enhancer region, contained matching profiles for 8-mers predicted from ldwgbm and npwgbm that are gc rich and similar to the gbm k-mers however, mutation of the single gli tfbs ablates the hh response, indicating that the presence of this gbm is required to transduce hh signaling. mutation of the gbms annotated by the green boxes for hhip, hipk ptch and scube is also sufficient to abrogate hh signal transduction. for boc, dpp and tgfbr which show enhancer activity that remains after ablation of the gbm, there were no sequence characteristics that were indicative of a shared feature responsible for this remaining response to induction by gli the boc profile was the only one that contained a high weighted k-mer that was unique to the ldwgbm profile. this k-mer was similar to a krox motif and may be enriched in the ldwgbm dataset as a context specific transcription factor, since it has roles in limb development. in general, most of the tested regions contain distributed high weighted 8-mers in addition to the central gbm and had profiles that contained consistent peaks in both ldwgm and npwgbm datasets.fig. k-mer weights plotted across sequences that show enhancer activity. diagrams were generated in ucsc genome browser and show coordinate information for regions annotated to ptch, boc, dpp, gli, hhip, hipk, ptch, scube, shh and tgfbr. green boxes represent gbm. weights for ldwgbm and npwgbm are represented by the red and yellow lines, respectively. refseq gene annotations are represented in blue. a putative krox- tfbs that has a high weight in the ldwgbm classifier but not the npwgbm classifier occurs in the sequence annotated to boc. note that most sequences show weighted k-mers located several hundred bp from the central gbm, suggesting that sequence motifs that predict hh enhancer activity may be distributed  we have utilized the kmer-svm machine learning approach to examine four existing gli chip databases and to generate classifiers that can empirically predict functional hh enhancers from genomic dna. the analysis was facilitated by a new gbs library generated from a meta-analysis of genomic regions identified in in vivo binding studies. when compared to the previous library derived from in vitro binding studies, this new gbs library resulted in the identification of a subset of potential gli transcription factor binding sites across the mouse genome. of nine predicted hh target regulatory regions tested, seven were functionally verified as gli-dependent. all of the tested regions were annotated to genes believed to be components of the hh pathway and important determinants of the hh response. given the high success rate of hh enhancer prediction in this small scale screen, it is quite possible that a large number of the other predicted regions may harbor hh enhancer activity. 
however, surprisingly few of the critical enhancer elements that control these pathway target genes have been identified despite the fact that such elements are central determinants of hh signaling activity. using these datasets, we tested whether a meta-analysis of gli binding sites, coupled with a machine learning approach, could reveal genomic features that could be used to empirically identify hh-regulated enhancers linked to loci of the hh signaling pathway. the kmer-svm machine learning approach identifies those features and can successfully predict functional hh regulatory regions in genomic dna surrounding hh pathway molecules and likely, other hh targets. recently, chip studies have been carried out in multiple tissue contexts using mouse models carrying flag-tagged gli proteins. additionally, the library of enriched gli binding motifs that we have identified may allow improved identification of functional gli binding sites. the results suggest that hh enhancer regions share common sequence features. for functional analysis, we selected nine regions which were annotated to putative hh pathway molecules and found that seven exhibited gli-dependent activity, indicating that they are directly regulated by hh signaling. importantly, a number of genes that encode hh pathway molecules are themselves targets of hh signaling, allowing pathway regulation by an intricate balance of feed-back activation and inhibition. a meta-analysis of four existing gliflag datasets revealed a library of gli binding motifs that was substantially more restricted than the potential sites predicted by previous in vitro binding studies. these enhancers integrate the hh signal in a context-specific manner to control the spatiotemporal pattern of target gene expression. the hedgehog signaling pathway, acting through three homologous transcription factors in vertebrates, plays multiple roles in embryonic organ development and adult tissue homeostasis. a machine learning method was then applied to these datasets and enriched k-mers were identified that, when applied to the mouse genome, predicted as many as potential hh enhancers.
 the hedgehog signaling pathway, acting through three homologous transcription factors in vertebrates, plays multiple roles in embryonic organ development and adult tissue homeostasis. at the level of the genome, gli factors bind to specific motifs in enhancers, some of which are hundreds of kilobases removed from the gene promoter. these enhancers integrate the hh signal in a context-specific manner to control the spatiotemporal pattern of target gene expression. importantly, a number of genes that encode hh pathway molecules are themselves targets of hh signaling, allowing pathway regulation by an intricate balance of feed-back activation and inhibition. however, surprisingly few of the critical enhancer elements that control these pathway target genes have been identified despite the fact that such elements are central determinants of hh signaling activity. recently, chip studies have been carried out in multiple tissue contexts using mouse models carrying flag-tagged gli proteins. using these datasets, we tested whether a meta-analysis of gli binding sites, coupled with a machine learning approach, could reveal genomic features that could be used to empirically identify hh-regulated enhancers linked to loci of the hh signaling pathway.  a meta-analysis of four existing gliflag datasets revealed a library of gli binding motifs that was substantially more restricted than the potential sites predicted by previous in vitro binding studies. a machine learning method was then applied to these datasets and enriched k-mers were identified that, when applied to the mouse genome, predicted as many as potential hh enhancers. for functional analysis, we selected nine regions which were annotated to putative hh pathway molecules and found that seven exhibited gli-dependent activity, indicating that they are directly regulated by hh signaling.  the results suggest that hh enhancer regions share common sequence features. the kmer-svm machine learning approach identifies those features and can successfully predict functional hh regulatory regions in genomic dna surrounding hh pathway molecules and likely, other hh targets. additionally, the library of enriched gli binding motifs that we have identified may allow improved identification of functional gli binding sites. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords hedgehog signalingenhancersmachine learninggliuniversity of michigan nih training program in basic and translational digestive sciencet dk094775gurdziel katherine niddk nih hhsp dk062041gumucio deborah l. issue-copyright-statement© the author 2016 