BACKGROUND
regulation of transcription in eukaryote cells is controlled by the binding of transcription factors  to specific binding sites in the regulatory regions of their target genes.
in this way, transcription factor binding sites  play an important role in the regulation of gene expression.
unfortunately, tfbss are hard to identify; experimental approaches are laborious and costly, and computational analyses are plagued by high false positive rates.
tfbss are extremely short compared to the regions in which they are present, and tfs typically bind to a variety of motifs.
one of the many difficulties faced by tfbs detection approaches is that some tfbss are restricted in their location with regard to the transcription start site .
an extreme example of such a motif is the tata-box, which is positioned about  <dig> to  <dig> bases upstream of the tss.
however, computational analyses usually use sequences of a fixed length .
in such cases, the region in which genuine regulatory motifs are positioned is small compared to the input sequence length, which makes position-restricted tfbss hard to detect using standard approaches.
a number of studies have investigated the positional tendencies of nucleotide motifs and tfbss in genome-wide or large sets  of promoter sequences  <cit> .
approaches aiming at predicting locally enriched tfbss in smaller sets of sequences have also been reported.
many of these approaches involve counting the number of tfbss in windows within the regulatory regions, either using a sliding window approach  <cit>  or by binning the tfbss according to their position  <cit> .
one weakness of these approaches is that they process the tfbss in a binary way; either a site is present within the window, or it is not.
within the region defined by the window or bin itself the distribution of sites is ignored, and sites at the edge of the region contribute as much to the score as do sites located at the center.
this causes such approaches to be highly dependent on window sizes, and – especially for smaller window sizes – susceptible to noise: moving the region or changing its size even just a few bases can drastically change the score.
alternative methods for combining positional information and nucleotide motifs, such as one based on positional word counting  <cit> , are not suitable for application on smaller datasets where counts of tfbss tend to be too low to process in such a way.
here we introduce a novel method for the prediction of locally enriched tfbss in sets of promoter sequences of co-regulated genes.
our approach is based on the parzen window technique for density estimation, and is capable of converting predicted tfbss in even small sets of promoters into continuous functions that reflect local tfbs enrichment.
usage of different window function widths allows for the detection of both tfbs motifs with loosely defined positional preferences as well as tfbss that require very precise positioning.
furthermore, a random sampling technique is used to incorporate genome-wide tendencies into the evaluation of significance of local enrichment.
application on artificial and real datasets showed that our method compares favourably to other methods and measures.
we also applied our method to a number of promoter sequence sets regulating tissue- or condition-specific expression profiles, and detected regions with local enrichment of tfbss that could not be detected using a more standard tfbs over-representation analysis.
we illustrate the validity of our findings and how our approach can be used for further dissection of the architecture of regulatory regions.
methods
our proposed method for detection of locally enriched regulatory motifs is described below.
a workflow of our approach as well as details not covered below can be found in the supplementary methods section and figure s <dig> .
parzen window-based scoring of local tfbs enrichment
we start with a set s of n promoter sequences, indexed by s =  <dig>   <dig>  …, n. for each sequence s its ms associated tfbss are represented by their location relative to a landmark within each sequence,
ls=ls, <dig> ls, <dig> …,ls,ms.
in this study the landmark is the tss associated with each promoter sequence, and our goal is to find a region relative to the tss within these promoter sequences that is significantly enriched in tfbss.
we start the description of our approach by temporarily focusing on a window of size h, which is used to scan the set of sequences.
the probability p that a tfbs will fall within a region r can be expressed as:
p=∫rpxdx
the density function p can be estimated empirically by
px=1n∑s=1n∑i=1msϕx,ls,i,h
where n is the total number of tfbss in the set of n sequences.
here ϕ is a kernel function referred to as the parzen window function.
in a simple sliding window approach ϕ could be defined as follows:
ϕx,l,h={1hifx-l≤h/20otherwise
in this case p simply reflects the count of sites in the window around x.
in our approach, a first adjustment we made is to replace the simple parzen window function of eq.
<dig> with a smoother window function:
ϕx,l,h=12πh/22exp-l-x22h/ <dig>
that is, ϕ is a gaussian function fit over each position x in the set of promoters, with the distance to predicted tfbss determining their contributions to the score of x. the h value here is essentially a measure of the width of this gaussian.
other window functions can easily be introduced.
secondly, we adjusted the parzen window approach so that p is no longer an estimate of the density of tfbss, but of the local enrichment of tfbss in the region around position x in the promoter sequences.
we define ps to be the contribution of the tfbss of sequence s to the local enrichment score at position x;
psx=∑i=1msϕx,ls,i,h
and the total local enrichment score  at position x over the entire set of sequences s is
slocalx=∑s=1n1zspsx
here, zs is a normalization factor inspired by the zero or one occurrence per sequence  model that is frequently used in motif detection programs  <cit> :
zs=max <dig> ∑x=xstartxstoppsx
where xstart and xstop indicate the start and stop coordinates of the region of interest relative to the tss, respectively.
in practice, if a sequence contains multiple predicted tfbss, zs will be greater than  <dig>  resulting in a reduced contribution per site in this sequence to the slocal.
this scaling factor was introduced to limit the effect of a single sequence containing high numbers of predicted tfbss, especially in simple repeat sequences .
note that the slocal is not a probability function and does not sum to unity.
using equations  <dig> to  <dig>  a local enrichment score can be calculated for each tfbs motif at each bp in the region from xstart to xstop.
in practice we focused on the region from -2 kb to 1 kb around the tss, but to avoid irregularities at the boundaries of this region we based calculations on tfbss predicted in the region from -3 kb to +2 kb.
wider window functions lead to broad, slowly changing slocal values which are useful for detecting over-representation of tfbss with relatively loosely defined regions of preference.
narrower window functions on the other hand result in sharper peaks in slocal values, allowing for the detection of local enrichment in narrow, precisely defined regions.
in order to detect both broad and narrow regions of enrichment we used the values h/ <dig> =  <dig>   <dig>   <dig>   <dig>  and 200 bps.
evaluation of significance of local enrichment
the significance of observed slocal values was evaluated using a random sampling approach, in which we sample n sequences with their predicted binding sites from the genomic set of promoters.
sampling was done such that sampled sequences had a similar gc content composition as the input sequences, in an effort to limit biases .
using equations  <dig> to  <dig>  for each position x, we calculated
slocalsampled.
we repeated this sampling a large number of times and for each position x we calculated two p-values: a position-dependent one, pdep, and a position-independent one, pind.
pdep is defined as the proportion of sampled sets where
slocalsampled ≥ slocal, and pind is defined as the proportion of sampled sets where max x') ≥ slocal.
the use of the combination of both thresholds ensures that enriched regions are enriched in s in comparison to the same region in sampled sets , but also that the regions have a certain degree of enrichment irrespective of their position to the tss .
in order to reduce false positive predictions caused by multiple testing, we defined the thresholds for pdep as a function of the window function width, with values of  <dig> , 5e- <dig>  2e- <dig>  1e- <dig>  and 5e- <dig> for scores calculated with h/2 =  <dig>   <dig>   <dig>   <dig>  and 10 bps respectively.
here, more stringent thresholds for the smaller peak widths reflect the increasing number of independent tests performed as peak widths decrease.
the pind threshold was set uniformly to  <dig> .
for every region within the sequences with start position x <dig> and end position x <dig> passing both p value restrictions we defined the corresponding region of local enrichment to be the region from x1 - h/ <dig> to x2 + h/ <dig>
removal of redundant locally enriched regions
pairs of overlapping enriched regions detected for the same motif with the same window function width were merged together to form one single region.
redundancy between enriched regions detected using different window function width values was removed such that for each set of overlapping regions only  <dig> representative region was retained.
for each region, the number of tfbss contained was counted.
next, for each set of overlapping regions we retained the largest region, unless there were smaller regions containing more than 2/ <dig> of the predicted tfbss in the largest region.
in the latter case we retained the smallest region containing more than 2/ <dig> of the predicted tfbss.
this approach ensured that, where possible, narrower enriched regions containing high numbers of tfbss were returned.
a random sampling approach that takes into account gc content biases
at several points in this study we evaluate the significance of findings using a random sampling strategy.
observations made for the set of sequences s were compared to values obtained in a large number of sets of n sequences randomly sampled from the genomic set of sequences in a way that limits the difference in gc content between sampled sequences and the sequences in s. the p-value of observed values was estimated by the ratio of sampled sets having a higher value than the observed one.
first, we clustered the genomic set of promoter sequences by their gc content.
each sequence was represented by a vector of  <dig> values representing the gc content in bins of 100 bps in the region from -1 kb to +1 kb.
values for each bin were scaled to have mean  <dig> and standard deviation  <dig>  clustering was done using k-means clustering using the hartigan-wong algorithm, with k =  <dig> to  <dig>  for each value of k, clustering was done  <dig> times using random initializations, each run had up to  <dig> iterations, and the result with the smallest sum of squares between samples and assigned centers was retained.
the clusters obtained with k =  <dig> corresponded to sequences with high  and low gc content, respectively.
clusters obtained by clustering with k >  <dig> corresponded to variations of both classes .
in our sampling approach, given a set s of n input sequences and a value of k*, we get the number of sequences in s belonging to each cluster,
c <dig> c <dig> …,ck*, with
∑i=1k* ci = n, and randomly sample ci sequences from each cluster of the genomic set of sequences.
here, k* represents the optimal k value, chosen for each set of input sequences s in a way that limited the difference in gc content profiles between s and sampled sequences.
for determining k*, for each k value , we randomly sampled ci sequences from each cluster of the genomic set of sequences, and, for the thus obtained n sampled sequences, calculated the average gc content in each bin of 100 bps.
the root-mean-square deviation  between these gc content values and the average gc content values of the sequences in s is calculated.
this sampling was repeated  <dig> times, allowing us to calculate the average rmsd and its standard deviation  .
the estimate of k* was then:
k*=argminkk|rmsdk≤rmsdk+1+sdk+ <dig>
which is inspired by the gap statistic proposed to estimate the optimal k in k-means clustering  <cit> .
over-representation index
for the purpose of comparison, as an indication of global enrichment  of a tfbs in s we employed the over-representation index   <cit> , based on tfbss predicted in the 1000 bps upstream of tsss.
ori=densitysdensitygenomic×proportionsproportiongenomic
where densitys is the number of predicted tfbss in s divided by n, densitygenomic the number of predicted tfbss in the genomic set of promoter sequences divided by the total number of promoter sequences, and proportions and proportiongenomic the ratio of sequences that have at least  <dig> predicted tfbs, in s and the genome-wide set of promoters, respectively.
the significance of observed ori values is evaluated using a sampling approach equivalent to the one used for local enrichment scores, with pori being the proportion of random samples with higher ori scores than the value observed in s. see  for a flowchart of this approach.
after application on artificial datasets, and based on our experience in previous analyses, we set the threshold for pori to  <dig> .
comparison to existing methods and measures
we compared the performance of our approach to that of existing methods and measures for the prediction of local and global enrichment, using artificial as well as real data.
a number of realistic, artificial sequence datasets were prepared in which instances of known regulatory motifs were inserted into specific regions with respect to the tsss.
real datasets consisted of  <dig> sets of genes with similar expression profiles in dendritic cells upon stimulation with lps, for which chip-seq time-course data is available for  <dig> tfs  <cit> .
using the chip-seq data we inferred the tfs controlling the expression of each set of genes.
on these datasets we applied the methods and evaluated their performance in terms of precision, recall, and f-measure.
for a more detailed description of the construction of the artificial and real datasets, the methods and measures we included in the comparison, and the measures of performance we used, we refer to the supplementary methods section .
summary of additional materials and methods
for a more detailed description we refer to the supplementary methods section .
in brief, sets of co-expressed genes were defined based on micro-array gene expression data for  <dig> human and  <dig> mouse tissues and cell types from the gnf geneatlas dataset  <cit> , and for mouse dendritic cells  stimulated with a number of immune stimuli  <cit> .
for all genes of the human  and mouse  genome, we scanned regions from -3 kb to +2 kb around the tss for sites for a set of non-redundant, vertebrate position weight matrices  constructed from transfac  <cit>  and jaspar  <cit>  pwms.
analysis of local enrichment of tfbss was based on these predicted sites.
evolutionary conservation of tfbss within enriched regions was evaluated using phastcons scores as available on the ucsc genome browser  <cit> .
z scores for phastcons scores corresponding to bases included in tfbss were calculated based on the average and standard deviation of phastcons scores of an equal number of bases located at equal distances to randomly selected tsss.
enrichment of weak tfbss in locally enriched regions was evaluated by comparing the count of weak sites within predicted regions in s with the count in randomly sampled sets of sequences.
RESULTS
parzen window-based detection of local enrichment of tfbss
the starting point for many tfbs analyses is a set of predicted tfbs in a set of regulatory sequences believed to be under the control of the same regulatory mechanism.
finding local over-representation of tfbss in a set of sequences can be thought of as finding a region at a certain distance from a landmark, such as the tss, in which there is a significant enrichment or increase in density of tfbss compared to a reference set of sequences.
various general approaches for the estimation of densities have been described, such as frequency histograms and kn-nearest-neighbor estimates  <cit> .
here, we focused our attention on the parzen window approach, a non-parametric technique where the density estimate is obtained by summing over all samples weighted by a function of choice, the so-called window function  <cit> .
parzen widow approaches have been used in bioinformatics for chip-seq peak calling  <cit> , but to our best knowledge it has never been used for the analysis of regulatory motifs.
in our analysis, samples consist of predicted tfbss in promoter sequences of co-expressed genes, and as a window function we use a gaussian function of the distance to each tfbs.
we choose a gaussian function because it results in a smooth enrichment score, and because it has only one additional parameter to specify: the width of the gaussian.
we made a number of adjustments to the original parzen window approach , which allow us to estimate local enrichment scores  for each tfbs.
a sampling approach is then used to estimate the significance of observed scores in a way that takes into account gc content profiles of sequences.
to illustrate our approach we first focus our attention on the predicted sites for a selection of  <dig> pwms in the promoters of  <dig> mouse genes with tissue-specific expression in small and large intestine.
a first important point is that for each of these  <dig> motifs, a region of local enrichment was found , while only one of them  was found to be significantly over-represented by a standard over-representation approach.
secondly, between the  <dig> motifs there is a diversity in the width of the region of enrichment : the tata box and rxr/rar/vdr motifs show very narrow slocal peaks, while the curves of the hnf <dig> and especially the bhlh motifs  show very broad peaks.
for comparison, the curves for expected slocal values based on the genome-wide set of promoter sequences are shown in figure 2b.
the peak for the tata box corresponds to its known location around positions - <dig> to - <dig>  other motifs show a smaller degree of variation in function of the position relative to the tss, and some have a more or less uniform distribution of sites.
in order to take into account such genome-wide tendencies, we employ a sampling approach for the evaluation of significant local enrichment .
figure 2c shows the position dependent p-values  estimated by this sampling approach.
importantly, the fact that the tata box is found to be significantly enriched in the promoters of intestine-specific genes implies that this motif is significantly enriched in the region just upstream of the tss, even taking into account its strong genome-wide preference for this region , suggesting that the tata box plays a role in the regulation of these genes.
in this way, our approach is able to find regulatory motifs that are missed by standard approaches, it can detect both broad and narrow regions of local enrichment of motifs, and it takes into account genome-wide tendencies during the evaluation of significance of enrichment.
several studies have indicated differences in gc content and cpg scores between housekeeping genes and tissue- or condition-specific genes  <cit> , and a recent study has suggested that analysis of regulatory regions can lead to better results when treating cpg island-associated sequences and non-cpg island-associated sequences separately  <cit> .
since pwms with high  gc content naturally tend to have more hits in sequences with high  gc content, local regions of higher  gc content in an input set of sequences can easily result in apparent local enrichment of gc-rich  tfbss.
to avoid such biases we designed a way to evaluate the significance of peaks using a sampling approach in which sampled sequences have a gc content profile similar to that of the input sequences .
figure  <dig> illustrates the importance of this approach using the set of mouse promoters with high activity in small and large intestine as an example.
figure 1a shows the gc content profile in the intestine-specific promoters, as well as in sequences randomly sampled using k =  <dig>   <dig>   <dig> cluster.
if we were to ignore gc content profiles of sequences, sampling from the genomic set of sequences  would result in sampled sets having a significantly higher gc content than the input sequences, especially around the tss.
this in turn would lead to high numbers of false positive enriched regions for at-rich tfbss.
however, when clustering genomic sequences into  <dig> clusters according to their gc content profile , and sampling according to the number of input sequences in each cluster, the difference in average gc content between sampled and input sequences is strongly reduced.
doing the same using  <dig> clusters  does not further decrease the difference in this particular case.
figure 1b shows the discrepancy between gc content of the same set of intestine-specific promoters and sampled sets of promoters as a function of the number of clusters k used to sample.
although the rmsd has a tendency to decrease as k increases, higher k values also result in smaller sets of promoter sequences to sample from.
as a compromise, we therefore used the clusters obtained using k*, the smallest k value where the rmsd is smaller than the rmsd + sd of k +  <dig> clusters.
in the case of promoters with high expression in small and large intestine, the random sampling for the evaluation of significance of local enrichment of tfbss was thus done according to k* =  <dig>  since the existence of high-cg genes and low-cg genes is well established, k* =  <dig> was not considered.
the k* values were  <dig> in  <dig> out of  <dig> sets of co-expressed genes in mouse, and in  <dig> out of  <dig> sets in human, indicating that in many cases a simple distinction between gc-rich and gc-poor sequences is sufficient.
for a number of sets, however, high k* values were found .
comparison to existing methods and measures
we compared the accuracy of our method to that of other methods and measures for the prediction of local and global motif enrichment, using artificial sequence datasets and real datasets based on rna-seq and chip-seq data  <cit> .
both artificial and real datasets were constructed in a way that reflects a typical application of our method, e.g.
promoters of sets of genes thought to be under the control of similar regulatory mechanisms .
an overview of the precision, recall, and f-measure of each method and measure is shown in table  <dig>  the comparison shows that our method compares favourably to other methods.
recall is relatively high , while precision is the highest of the investigated methods for both artificial and real datasets .
interestingly, our method also had relatively high performance when using a simple uniform window function, though both recall and precision were lower than when using gaussian-based window function.
although a number of approaches had higher recall values than our method, these methods tend to also have considerably lower precision.
when methods for the prediction of global enrichment are run with more stringent thresholds in order to improve precision, recall drops and especially artificial tfbss inserted with lower insertion rates  tend to be missed .
other methods for the prediction of positional preferences or local enrichment of tfbss suffer in general from low precision, although we have to point out that several methods and measures were originally designed for slightly different purposes than the one investigated here.
the prediction of de novo motifs is in general regarded as a much harder problem than the detection of enrichment of known motifs.
this is reflected in the low recall and low precision of a-glam and free  <cit> , two methods which do not use an input set of pwms, but are based on the analysis of k-mers.
for each method or measure the type of measure , the recall, precision, and f-measure is given for the artificial and real datasets, as well as a reference.
methods are sorted by decreasing f-measure obtained on the real datasets.
p value threshold  <dig> ;  p value threshold  <dig> ;  p value threshold 1e-4;  sig threshold 0;  sig threshold  <dig>
in addition, we observed how our method significantly increases the accuracy of pwm-based tfbs predictions.
when we assume that only predicted tfbs present inside bound regions of the tf in question  are truly functional, and regard any other predicted tfbss as false positives, pwm-based tfbs predictions resulted in a recall of  <dig> %, with a low precision of  <dig> %.
however, when we filtered the pwm-based predicted tfbss using the predicted regions of enrichment for the regulatory motif in question, while recall dropped  <dig> -fold to  <dig> %, there was a more than 7-fold increase in precision to  <dig> %, resulting from a strongly reduced number of false positive predictions.
application on a large number of gene sets
next, we applied our approach on a large number of sets of genes with tissue- and cell type-specific expression, in human and mouse  <cit> .
genes were clustered into  <dig> and  <dig> clusters of co-expressed genes for the mouse and human case, respectively.
promoter regions for all genes were scanned using a set of  <dig> non-redundant vertebrate pwms, and we predicted local tfbs motif enrichment in each of the clusters using our approach, as well as tfbs enrichment using a standard motif over-representation approach.
in total our approach predicted  <dig> and  <dig> regions with local enrichment of a regulatory motif in the mouse and human gnf geneatlas datasets, respectively .
the fact that our approach could find locally enriched motifs that could also be detected using standard approaches demonstrates the robustness of our approach .
for example, sites for two ets family tfs, including pu.
<dig>  and for nf-κb, and irfs were found in various immune-related cell types  <cit> , and sites for hnf <dig> and hnf <dig> in liver-specific and liver/kidney-specific genes  <cit> .
moreover, table  <dig> shows a small selection of locally enriched motifs that could not be detected by the standard motif over-representation approach for the mouse clusters .
for example, we found enrichment of binding sites of ets family tfs in promoters driving high expression in mouse b cells and t cells, and a peak for a creb binding motif in a set of promoters associated with genes with preferential expression in mouse testis.
note that in many cases, motifs with local enrichment also show a certain degree of general enrichment , but on a level that cannot be regarded as significant .
these results clearly indicate that our approach is able to detect regulatory features that are missed by a standard motif over-representation analysis, and that in many cases the detected regulatory motifs are known to play a role in the dataset in which they were detected.
a selection of regulatory motifs is shown for which regions of local enrichment were detected in mouse tissues and cell types of the gnf geneatlas dataset.
the tissue, the start and stop position of the region, the h/ <dig> used, the ori p-value, and references supporting the role of the regulatory motif in the tissue in question are shown.
general tendencies of local tfbs enrichment
in general, locally enriched regions were found to be present proximal to the tss, roughly in the region - <dig> to + <dig> with a peak around position - <dig> .
interestingly, in the human case we observed a region that is relatively poor in enriched regions, roughly between positions - <dig> to - <dig>  suggesting an upper limit for distances over which positioning of tfbs relative to the tss is of biological importance.
in mouse too, the number of times a nucleotide was included in an enriched region dropped to a low basal level for positions upstream of position - <dig>  a similar limit was observed for the region downstream of the tss, although both in human and in mouse the slope of the curve is not as steep.
while the peak in the region just upstream of the tss seems highly significant and suggests a recommended region for regulatory motif analysis roughly between positions - <dig> and + <dig>  it is possible that it is partly caused by a bias of past studies to focus mainly on this region for the identification of tfbss.
this might in turn bias the pwm data we used for detection of local enrichment of tfbss.
on the other hand, given the flexibility of the dna double helix it is unlikely that very precise positioning of tfbss is necessary at larger distances from the tss.
regulatory motifs show a preference for similar regions across different sets of promoters
an additional indication of the biological validity of our findings is the fact that for many motifs we found local enrichment in similar regions in different clusters of genes.
one example is the tata box, which is known to be present at a specific distance from the tss and which we found to be significantly enriched in the region around positions - <dig> to - <dig> in two sets of promoters.
we also observed similar tendencies for motifs that are not generally known to be under strong positional restrictions.
some examples are ets family tfbss, which we detected concentrated around the region immediately upstream of the tss in  <dig> human and  <dig> mouse clusters , and nrf- <dig> binding sites, which we found in  <dig> human and  <dig> mouse clusters, here too in all cases in the region immediately upstream of the tss.
importantly, not all motifs were enriched in the upstream regions: zfx motifs, for example, were found to be enriched in regions roughly between  <dig> to  <dig> bases downstream of the tss in  <dig> clusters.
standard motif over-representation analyses are likely to miss such regulatory patterns, as this region is usually not included in the input.
region with local enrichment of ets domain tf binding motifs were found in different sets of sequences.
the species, and the tissues and cell types associated with the promoters in which the motif was found are listed, as well as the regions of local enrichment.
tfbss in locally enriched regions tend to have high evolutionary conservation
another illustration of the validity of our findings is the relatively high evolutionary conservation of tfbss present within enriched regions, compared to nucleotides within predicted tfbss at the same position relative to tsss of randomly sampled genes .
among the  <dig> enriched regions containing ≥  <dig> tfbss that were detected in the human gnf geneatlas sets,  <dig>  contain tfbss that are on average more conserved  than randomly picked positions, while only  <dig>  contain sites that are on average less conserved .
similarly, for mouse datasets, we found  <dig> out of  <dig>  regions with highly conserved sites, and only  <dig>  with low conservation .
in both human and mouse, the vast majority of the peaks with high conservation were located proximal to the tss .
among the relatively more conserved sites are the predicted sites for yy <dig> in genes specifically expressed in mouse thymus.
the bases corresponding to the  <dig> yy <dig> sites predicted in the region - <dig> to + <dig> have an average phastcons score of  <dig>  .
evolutionary conservation was not only found on the level of predicted tfbss, but also on the level of predicted regions of local enrichment, as illustrated in figure  <dig> for promoters of testis-specific genes.
for the sites of pwms representing rfx <dig>  rfx tfs in general, a creb motif, and myb family tfs, similar regions of enrichment were predicted in human and mouse sequences.
locally enriched regions tend to be enriched for weak tfbss
recently a number of studies, including computational studies using thermodynamic models, have illustrated the importance of weak tfbss in the regulation of transcription  <cit> .
although pwm scores are believed to be correlated with the tf-dna binding energy  <cit> , the sheer number of sites having only moderately high pwm scores prohibits most computational methods from analyzing them in more detail.
one advantage our approach offers is that we can restrict the region of interest to the region preferred by the stronger tfbss, and evaluate if weaker binding sites show a similar local enrichment .
our analysis of weak tfbss in the enriched regions found in the human clusters shows that  <dig> out of  <dig> regions  show a significant enrichment  of weaker tfbss .
similarly, in mouse datasets,  <dig> out of  <dig> regions  were enriched for weak tfbss .
this result does not include the presence of stronger tfbss, nor weaker sites overlapping with them.
for example, promoters of genes expressed in human liver are enriched for hnf <dig> tfbss in the region - <dig> to + <dig>  our analysis indicates that in addition to the  <dig> strong hnf <dig> sites, there are  <dig> weak tfbss predicted in the region of local enrichment in the promoters in this set of genes, .
application on sets of genes with differential expression upon tlr stimulation
the regulatory pathways controlling toll-like receptor  signalling have been relatively well studied, and main regulators such as nf-κb and irfs have been identified.
however, it is likely that promoters of genes that are induced or repressed after tlr stimulation contain additional regulatory motifs that are still to be identified.
from gene expression data taken from mouse dcs stimulated with  <dig> tlr ligands  <cit> , we prepared  <dig> clusters of co-expressed genes.
applying our method to these clusters, we could detect locally enriched regions for the sites of creb-binding tfs, nf-κb, and irfs in clusters of transiently induced genes with peak expression around  <dig>   <dig>  and 6 hrs after stimulation, respectively .
these motifs could also be found by standard motif over-representation analysis.
however, our approach also detected a local enrichment for stat binding sites in a cluster with induced expression peaking around 6 hours after stimulation  <cit> , which could not be detected by the standard approach.
in addition, we could make a number of findings that might offer hypotheses for further wet-lab experiments.
one is the enrichment of tata box motifs only in the set of promoters corresponding to genes with the fastest transient activation , which suggests that the tata box plays a role in the rapid induction of these genes.
finally, a number of enriched regions downstream of the tss were predicted.
these include regions enriched in sites for hmx <dig> , rfx <dig> sites , mtf <dig> , zfp <dig> and e-box motifs .
CONCLUSIONS
although various studies have illustrated that in eukaryotic genomes cis-regulatory motifs can be positioned several kilobases or even megabases away from their target genes, it has also been reported that in a number of cases tfbss show a tendency to be present at a more or less fixed distance with regard to the tss.
nevertheless, in general, no clear positional preferences have been described for most regulatory motifs, even though such information could be extremely useful for their prediction.
here, we present a novel method for detecting locally enriched tfbss in the regulatory regions of sets of co-regulated genes.
our approach is based on a non-parametric approach for sample density estimation, with adjustments which allow it to detect sequence regions that have a significant local enrichment in tfbss, on a 1 bp resolution.
an online tool of our approach, which we call locamo finder , is available at http://sysimm.ifrec.osaka-u.ac.jp/tfbs/locamo/.
our approach evaluates both positioning and enrichment of tfbss simultaneously, using a set of control sequences as a reference.
we implemented our method for the estimation of significance of enrichment in a way that takes into account gc content profiles of the input set of promoter sequences.
evaluation of significance is done against randomly selected promoters that have gc content profiles similar to the input sequences.
we showed how our approach has a high accuracy compared to other methods and measures for local and global motif enrichment.
we applied our method on a large number of sets of genes with tissue- or cell type-specific expression, as well as on a number of sets of genes with similar expression profiles after tlr ligand stimulation in mouse dcs.
detected regions of local enrichment of tfbss are supported by known regulatory interactions reported in literature, as well as by the observation that several regulatory motifs are found to be locally enriched in similar regions in different sets of promoters.
in addition, detected regions tend to contain tfbss with higher evolutionary conservation than expected, and they also tend to be enriched for weak tfbss.
together, these results illustrated the usefulness and validity of our approach.
approaches for finding local enrichment of tfbss heavily rely on tss annotations.
recent studies have shown that a large fraction of genes have several tsss, and that different tsss might allow different amounts of variation in the bases from which transcription is initiated  <cit> .
future approaches that incorporate such features are likely to present more biological insights into the relationship between tfbs positioning and transcription initiation and tss variability.
abbreviations
dc: dendritic cell; ori: over-representation index; pwm: position weight matrix; rmsd: root-mean-square deviation; sd: standard deviation; tf: transcription factor; tfbs: transcription factor binding site; tlr: toll-like receptor; tss: transcription start site; zoops: zero or one occurrence per sequence
competing interests
the authors declare that they have no competing interests.
authors’ contributions
av conceived of the study and performed the bioinformatics analysis, and prepared the manuscript.
yk, st, and kma contributed to statistical analysis, discussion and interpretation of results and with drafting the manuscript.
sa and dms supervised the project and helped with discussion and interpretation of results and with drafting the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
supplementary text and figures.
a file containing supplementary material, including detailed description of methods and results, as well as supplementary figures and tables.
click here for file
additional file 2
supplementary tables s <dig>  s <dig>  and s <dig>  a spreadsheet file containing predicted regions of local enrichment and globally enriched tfbs motifs for human and mouse gnf geneatlas datasets, and for the sets of promoters obtained from tlr-stimulated mouse dc cells.
click here for file
acknowledgements
the authors would like to thank all members of the laboratory of systems immunology and the laboratory of host defense for helpful discussions and advice; dr k nakai for helpful suggestions; m ogawa, e kamada and m kageyama for secretarial assistance.
we also thank the reviewers for their suggestions, which improved the paper considerably.
funding
this work was supported by the special coordination funds of the japanese ministry of education, culture, sports, science and technology, and the ministry of health, labour and welfare in japan, the japan society for the promotion of science through funding program for world-leading innovative r&d on science and technology , and by a kakenhi grant-in-aid for scientific research  from the japan society for the promotion of science.

 regulation of transcription in eukaryote cells is controlled by the binding of transcription factors to specific binding sites in the regulatory regions of their target genes. in this way, transcription factor binding sites play an important role in the regulation of gene expression. unfortunately, tfbss are hard to identify; experimental approaches are laborious and costly, and computational analyses are plagued by high false positive rates. tfbss are extremely short compared to the regions in which they are present, and tfs typically bind to a variety of motifs. one of the many difficulties faced by tfbs detection approaches is that some tfbss are restricted in their location with regard to the transcription start site. an extreme example of such a motif is the tata-box, which is positioned about to bases upstream of the tss. however, computational analyses usually use sequences of a fixed length. in such cases, the region in which genuine regulatory motifs are positioned is small compared to the input sequence length, which makes position-restricted tfbss hard to detect using standard approaches. a number of studies have investigated the positional tendencies of nucleotide motifs and tfbss in genome-wide or large sets of promoter sequences. approaches aiming at predicting locally enriched tfbss in smaller sets of sequences have also been reported. many of these approaches involve counting the number of tfbss in windows within the regulatory regions, either using a sliding window approach or by binning the tfbss according to their position. one weakness of these approaches is that they process the tfbss in a binary way; either a site is present within the window, or it is not. within the region defined by the window or bin itself the distribution of sites is ignored, and sites at the edge of the region contribute as much to the score as do sites located at the center. this causes such approaches to be highly dependent on window sizes, and – especially for smaller window sizes – susceptible to noise: moving the region or changing its size even just a few bases can drastically change the score. alternative methods for combining positional information and nucleotide motifs, such as one based on positional word counting, are not suitable for application on smaller datasets where counts of tfbss tend to be too low to process in such a way. here we introduce a novel method for the prediction of locally enriched tfbss in sets of promoter sequences of co-regulated genes. our approach is based on the parzen window technique for density estimation, and is capable of converting predicted tfbss in even small sets of promoters into continuous functions that reflect local tfbs enrichment. usage of different window function widths allows for the detection of both tfbs motifs with loosely defined positional preferences as well as tfbss that require very precise positioning. furthermore, a random sampling technique is used to incorporate genome-wide tendencies into the evaluation of significance of local enrichment. application on artificial and real datasets showed that our method compares favourably to other methods and measures. we also applied our method to a number of promoter sequence sets regulating tissue- or condition-specific expression profiles, and detected regions with local enrichment of tfbss that could not be detected using a more standard tfbs over-representation analysis. we illustrate the validity of our findings and how our approach can be used for further dissection of the architecture of regulatory regions. methods our proposed method for detection of locally enriched regulatory motifs is described below. a workflow of our approach as well as details not covered below can be found in the supplementary methods section and figure s. parzen window-based scoring of local tfbs enrichment we start with a set s of n promoter sequences, indexed by s =     …, n. for each sequence s its ms associated tfbss are represented by their location relative to a landmark within each sequence, ls=ls, ls, …,ls,ms. in this study the landmark is the tss associated with each promoter sequence, and our goal is to find a region relative to the tss within these promoter sequences that is significantly enriched in tfbss. we start the description of our approach by temporarily focusing on a window of size h, which is used to scan the set of sequences. the probability p that a tfbs will fall within a region r can be expressed as: p=∫rpxdx the density function p can be estimated empirically by px=1n∑s=1n∑i=1msϕx,ls,i,h where n is the total number of tfbss in the set of n sequences. here ϕ is a kernel function referred to as the parzen window function. in a simple sliding window approach ϕ could be defined as follows: ϕx,l,h={1hifx-l≤h/20otherwise in this case p simply reflects the count of sites in the window around x. in our approach, a first adjustment we made is to replace the simple parzen window function of eq. with a smoother window function: ϕx,l,h=12πh/22exp-l-x22h/ that is, ϕ is a gaussian function fit over each position x in the set of promoters, with the distance to predicted tfbss determining their contributions to the score of x. the h value here is essentially a measure of the width of this gaussian. other window functions can easily be introduced. secondly, we adjusted the parzen window approach so that p is no longer an estimate of the density of tfbss, but of the local enrichment of tfbss in the region around position x in the promoter sequences. we define ps to be the contribution of the tfbss of sequence s to the local enrichment score at position x; psx=∑i=1msϕx,ls,i,h and the total local enrichment score at position x over the entire set of sequences s is slocalx=∑s=1n1zspsx here, zs is a normalization factor inspired by the zero or one occurrence per sequence model that is frequently used in motif detection programs: zs=max ∑x=xstartxstoppsx where xstart and xstop indicate the start and stop coordinates of the region of interest relative to the tss, respectively. in practice, if a sequence contains multiple predicted tfbss, zs will be greater than resulting in a reduced contribution per site in this sequence to the slocal. this scaling factor was introduced to limit the effect of a single sequence containing high numbers of predicted tfbss, especially in simple repeat sequences. note that the slocal is not a probability function and does not sum to unity. using equations to a local enrichment score can be calculated for each tfbs motif at each bp in the region from xstart to xstop. in practice we focused on the region from -2 kb to 1 kb around the tss, but to avoid irregularities at the boundaries of this region we based calculations on tfbss predicted in the region from -3 kb to +2 kb. wider window functions lead to broad, slowly changing slocal values which are useful for detecting over-representation of tfbss with relatively loosely defined regions of preference. narrower window functions on the other hand result in sharper peaks in slocal values, allowing for the detection of local enrichment in narrow, precisely defined regions. in order to detect both broad and narrow regions of enrichment we used the values h/ = and 200 bps. evaluation of significance of local enrichment the significance of observed slocal values was evaluated using a random sampling approach, in which we sample n sequences with their predicted binding sites from the genomic set of promoters. sampling was done such that sampled sequences had a similar gc content composition as the input sequences, in an effort to limit biases. using equations to for each position x, we calculated slocalsampled. we repeated this sampling a large number of times and for each position x we calculated two p-values: a position-dependent one, pdep, and a position-independent one, pind. pdep is defined as the proportion of sampled sets where slocalsampled ≥ slocal, and pind is defined as the proportion of sampled sets where max x') ≥ slocal. the use of the combination of both thresholds ensures that enriched regions are enriched in s in comparison to the same region in sampled sets, but also that the regions have a certain degree of enrichment irrespective of their position to the tss. in order to reduce false positive predictions caused by multiple testing, we defined the thresholds for pdep as a function of the window function width, with values of, 5e- 2e- 1e- and 5e- for scores calculated with h/2 =  and 10 bps respectively. here, more stringent thresholds for the smaller peak widths reflect the increasing number of independent tests performed as peak widths decrease. the pind threshold was set uniformly to. for every region within the sequences with start position x and end position x passing both p value restrictions we defined the corresponding region of local enrichment to be the region from x1 - h/ to x2 + h/ removal of redundant locally enriched regions pairs of overlapping enriched regions detected for the same motif with the same window function width were merged together to form one single region. redundancy between enriched regions detected using different window function width values was removed such that for each set of overlapping regions only representative region was retained. for each region, the number of tfbss contained was counted. next, for each set of overlapping regions we retained the largest region, unless there were smaller regions containing more than 2/ of the predicted tfbss in the largest region. in the latter case we retained the smallest region containing more than 2/ of the predicted tfbss. this approach ensured that, where possible, narrower enriched regions containing high numbers of tfbss were returned. a random sampling approach that takes into account gc content biases at several points in this study we evaluate the significance of findings using a random sampling strategy. observations made for the set of sequences s were compared to values obtained in a large number of sets of n sequences randomly sampled from the genomic set of sequences in a way that limits the difference in gc content between sampled sequences and the sequences in s. the p-value of observed values was estimated by the ratio of sampled sets having a higher value than the observed one. first, we clustered the genomic set of promoter sequences by their gc content. each sequence was represented by a vector of values representing the gc content in bins of 100 bps in the region from -1 kb to +1 kb. values for each bin were scaled to have mean and standard deviation clustering was done using k-means clustering using the hartigan-wong algorithm, with k =  to for each value of k, clustering was done times using random initializations, each run had up to iterations, and the result with the smallest sum of squares between samples and assigned centers was retained. the clusters obtained with k =  corresponded to sequences with high and low gc content, respectively. clusters obtained by clustering with k >  corresponded to variations of both classes. in our sampling approach, given a set s of n input sequences and a value of k*, we get the number of sequences in s belonging to each cluster, c c …,ck*, with ∑i=1k* ci = n, and randomly sample ci sequences from each cluster of the genomic set of sequences. here, k* represents the optimal k value, chosen for each set of input sequences s in a way that limited the difference in gc content profiles between s and sampled sequences. for determining k*, for each k value, we randomly sampled ci sequences from each cluster of the genomic set of sequences, and, for the thus obtained n sampled sequences, calculated the average gc content in each bin of 100 bps. the root-mean-square deviation between these gc content values and the average gc content values of the sequences in s is calculated. this sampling was repeated times, allowing us to calculate the average rmsd and its standard deviation. the estimate of k* was then: k*=argminkk|rmsdk≤rmsdk+1+sdk+ which is inspired by the gap statistic proposed to estimate the optimal k in k-means clustering. over-representation index for the purpose of comparison, as an indication of global enrichment of a tfbs in s we employed the over-representation index, based on tfbss predicted in the 1000 bps upstream of tsss. ori=densitysdensitygenomic×proportionsproportiongenomic where densitys is the number of predicted tfbss in s divided by n, densitygenomic the number of predicted tfbss in the genomic set of promoter sequences divided by the total number of promoter sequences, and proportions and proportiongenomic the ratio of sequences that have at least predicted tfbs, in s and the genome-wide set of promoters, respectively. the significance of observed ori values is evaluated using a sampling approach equivalent to the one used for local enrichment scores, with pori being the proportion of random samples with higher ori scores than the value observed in s. see for a flowchart of this approach. after application on artificial datasets, and based on our experience in previous analyses, we set the threshold for pori to. comparison to existing methods and measures we compared the performance of our approach to that of existing methods and measures for the prediction of local and global enrichment, using artificial as well as real data. a number of realistic, artificial sequence datasets were prepared in which instances of known regulatory motifs were inserted into specific regions with respect to the tsss. real datasets consisted of sets of genes with similar expression profiles in dendritic cells upon stimulation with lps, for which chip-seq time-course data is available for tfs. using the chip-seq data we inferred the tfs controlling the expression of each set of genes. on these datasets we applied the methods and evaluated their performance in terms of precision, recall, and f-measure. for a more detailed description of the construction of the artificial and real datasets, the methods and measures we included in the comparison, and the measures of performance we used, we refer to the supplementary methods section. summary of additional materials and methods for a more detailed description we refer to the supplementary methods section. in brief, sets of co-expressed genes were defined based on micro-array gene expression data for human and mouse tissues and cell types from the gnf geneatlas dataset, and for mouse dendritic cells stimulated with a number of immune stimuli. for all genes of the human and mouse genome, we scanned regions from -3 kb to +2 kb around the tss for sites for a set of non-redundant, vertebrate position weight matrices constructed from transfac and jaspar pwms. analysis of local enrichment of tfbss was based on these predicted sites. evolutionary conservation of tfbss within enriched regions was evaluated using phastcons scores as available on the ucsc genome browser. z scores for phastcons scores corresponding to bases included in tfbss were calculated based on the average and standard deviation of phastcons scores of an equal number of bases located at equal distances to randomly selected tsss. enrichment of weak tfbss in locally enriched regions was evaluated by comparing the count of weak sites within predicted regions in s with the count in randomly sampled sets of sequences.  parzen window-based detection of local enrichment of tfbss the starting point for many tfbs analyses is a set of predicted tfbs in a set of regulatory sequences believed to be under the control of the same regulatory mechanism. finding local over-representation of tfbss in a set of sequences can be thought of as finding a region at a certain distance from a landmark, such as the tss, in which there is a significant enrichment or increase in density of tfbss compared to a reference set of sequences. various general approaches for the estimation of densities have been described, such as frequency histograms and kn-nearest-neighbor estimates. here, we focused our attention on the parzen window approach, a non-parametric technique where the density estimate is obtained by summing over all samples weighted by a function of choice, the so-called window function. parzen widow approaches have been used in bioinformatics for chip-seq peak calling, but to our best knowledge it has never been used for the analysis of regulatory motifs. in our analysis, samples consist of predicted tfbss in promoter sequences of co-expressed genes, and as a window function we use a gaussian function of the distance to each tfbs. we choose a gaussian function because it results in a smooth enrichment score, and because it has only one additional parameter to specify: the width of the gaussian. we made a number of adjustments to the original parzen window approach, which allow us to estimate local enrichment scores for each tfbs. a sampling approach is then used to estimate the significance of observed scores in a way that takes into account gc content profiles of sequences. to illustrate our approach we first focus our attention on the predicted sites for a selection of pwms in the promoters of mouse genes with tissue-specific expression in small and large intestine. a first important point is that for each of these motifs, a region of local enrichment was found, while only one of them was found to be significantly over-represented by a standard over-representation approach. secondly, between the motifs there is a diversity in the width of the region of enrichment: the tata box and rxr/rar/vdr motifs show very narrow slocal peaks, while the curves of the hnf and especially the bhlh motifs show very broad peaks. for comparison, the curves for expected slocal values based on the genome-wide set of promoter sequences are shown in figure 2b. the peak for the tata box corresponds to its known location around positions - to - other motifs show a smaller degree of variation in function of the position relative to the tss, and some have a more or less uniform distribution of sites. in order to take into account such genome-wide tendencies, we employ a sampling approach for the evaluation of significant local enrichment. figure 2c shows the position dependent p-values estimated by this sampling approach. importantly, the fact that the tata box is found to be significantly enriched in the promoters of intestine-specific genes implies that this motif is significantly enriched in the region just upstream of the tss, even taking into account its strong genome-wide preference for this region, suggesting that the tata box plays a role in the regulation of these genes. in this way, our approach is able to find regulatory motifs that are missed by standard approaches, it can detect both broad and narrow regions of local enrichment of motifs, and it takes into account genome-wide tendencies during the evaluation of significance of enrichment. several studies have indicated differences in gc content and cpg scores between housekeeping genes and tissue- or condition-specific genes, and a recent study has suggested that analysis of regulatory regions can lead to better results when treating cpg island-associated sequences and non-cpg island-associated sequences separately. since pwms with high gc content naturally tend to have more hits in sequences with high gc content, local regions of higher gc content in an input set of sequences can easily result in apparent local enrichment of gc-rich tfbss. to avoid such biases we designed a way to evaluate the significance of peaks using a sampling approach in which sampled sequences have a gc content profile similar to that of the input sequences. figure  illustrates the importance of this approach using the set of mouse promoters with high activity in small and large intestine as an example. figure 1a shows the gc content profile in the intestine-specific promoters, as well as in sequences randomly sampled using k =  cluster. if we were to ignore gc content profiles of sequences, sampling from the genomic set of sequences would result in sampled sets having a significantly higher gc content than the input sequences, especially around the tss. this in turn would lead to high numbers of false positive enriched regions for at-rich tfbss. however, when clustering genomic sequences into clusters according to their gc content profile, and sampling according to the number of input sequences in each cluster, the difference in average gc content between sampled and input sequences is strongly reduced. doing the same using clusters does not further decrease the difference in this particular case. figure 1b shows the discrepancy between gc content of the same set of intestine-specific promoters and sampled sets of promoters as a function of the number of clusters k used to sample. although the rmsd has a tendency to decrease as k increases, higher k values also result in smaller sets of promoter sequences to sample from. as a compromise, we therefore used the clusters obtained using k*, the smallest k value where the rmsd is smaller than the rmsd + sd of k +  clusters. in the case of promoters with high expression in small and large intestine, the random sampling for the evaluation of significance of local enrichment of tfbss was thus done according to k* =  since the existence of high-cg genes and low-cg genes is well established, k* =  was not considered. the k* values were in out of sets of co-expressed genes in mouse, and in out of sets in human, indicating that in many cases a simple distinction between gc-rich and gc-poor sequences is sufficient. for a number of sets, however, high k* values were found. comparison to existing methods and measures we compared the accuracy of our method to that of other methods and measures for the prediction of local and global motif enrichment, using artificial sequence datasets and real datasets based on rna-seq and chip-seq data. both artificial and real datasets were constructed in a way that reflects a typical application of our method, e.g. promoters of sets of genes thought to be under the control of similar regulatory mechanisms. an overview of the precision, recall, and f-measure of each method and measure is shown in table  the comparison shows that our method compares favourably to other methods. recall is relatively high, while precision is the highest of the investigated methods for both artificial and real datasets. interestingly, our method also had relatively high performance when using a simple uniform window function, though both recall and precision were lower than when using gaussian-based window function. although a number of approaches had higher recall values than our method, these methods tend to also have considerably lower precision. when methods for the prediction of global enrichment are run with more stringent thresholds in order to improve precision, recall drops and especially artificial tfbss inserted with lower insertion rates tend to be missed. other methods for the prediction of positional preferences or local enrichment of tfbss suffer in general from low precision, although we have to point out that several methods and measures were originally designed for slightly different purposes than the one investigated here. the prediction of de novo motifs is in general regarded as a much harder problem than the detection of enrichment of known motifs. this is reflected in the low recall and low precision of a-glam and free, two methods which do not use an input set of pwms, but are based on the analysis of k-mers. for each method or measure the type of measure, the recall, precision, and f-measure is given for the artificial and real datasets, as well as a reference. methods are sorted by decreasing f-measure obtained on the real datasets. p value threshold; p value threshold; p value threshold 1e-4; sig threshold 0; sig threshold in addition, we observed how our method significantly increases the accuracy of pwm-based tfbs predictions. when we assume that only predicted tfbs present inside bound regions of the tf in question are truly functional, and regard any other predicted tfbss as false positives, pwm-based tfbs predictions resulted in a recall of %, with a low precision of %. however, when we filtered the pwm-based predicted tfbss using the predicted regions of enrichment for the regulatory motif in question, while recall dropped -fold to %, there was a more than 7-fold increase in precision to %, resulting from a strongly reduced number of false positive predictions. application on a large number of gene sets next, we applied our approach on a large number of sets of genes with tissue- and cell type-specific expression, in human and mouse. genes were clustered into and clusters of co-expressed genes for the mouse and human case, respectively. promoter regions for all genes were scanned using a set of non-redundant vertebrate pwms, and we predicted local tfbs motif enrichment in each of the clusters using our approach, as well as tfbs enrichment using a standard motif over-representation approach. in total our approach predicted and regions with local enrichment of a regulatory motif in the mouse and human gnf geneatlas datasets, respectively. the fact that our approach could find locally enriched motifs that could also be detected using standard approaches demonstrates the robustness of our approach. for example, sites for two ets family tfs, including pu. and for nf-κb, and irfs were found in various immune-related cell types, and sites for hnf and hnf in liver-specific and liver/kidney-specific genes. moreover, table  shows a small selection of locally enriched motifs that could not be detected by the standard motif over-representation approach for the mouse clusters. for example, we found enrichment of binding sites of ets family tfs in promoters driving high expression in mouse b cells and t cells, and a peak for a creb binding motif in a set of promoters associated with genes with preferential expression in mouse testis. note that in many cases, motifs with local enrichment also show a certain degree of general enrichment, but on a level that cannot be regarded as significant. these results clearly indicate that our approach is able to detect regulatory features that are missed by a standard motif over-representation analysis, and that in many cases the detected regulatory motifs are known to play a role in the dataset in which they were detected. a selection of regulatory motifs is shown for which regions of local enrichment were detected in mouse tissues and cell types of the gnf geneatlas dataset. the tissue, the start and stop position of the region, the h/ used, the ori p-value, and references supporting the role of the regulatory motif in the tissue in question are shown. general tendencies of local tfbs enrichment in general, locally enriched regions were found to be present proximal to the tss, roughly in the region - to + with a peak around position -. interestingly, in the human case we observed a region that is relatively poor in enriched regions, roughly between positions - to - suggesting an upper limit for distances over which positioning of tfbs relative to the tss is of biological importance. in mouse too, the number of times a nucleotide was included in an enriched region dropped to a low basal level for positions upstream of position - a similar limit was observed for the region downstream of the tss, although both in human and in mouse the slope of the curve is not as steep. while the peak in the region just upstream of the tss seems highly significant and suggests a recommended region for regulatory motif analysis roughly between positions - and + it is possible that it is partly caused by a bias of past studies to focus mainly on this region for the identification of tfbss. this might in turn bias the pwm data we used for detection of local enrichment of tfbss. on the other hand, given the flexibility of the dna double helix it is unlikely that very precise positioning of tfbss is necessary at larger distances from the tss. regulatory motifs show a preference for similar regions across different sets of promoters an additional indication of the biological validity of our findings is the fact that for many motifs we found local enrichment in similar regions in different clusters of genes. one example is the tata box, which is known to be present at a specific distance from the tss and which we found to be significantly enriched in the region around positions - to - in two sets of promoters. we also observed similar tendencies for motifs that are not generally known to be under strong positional restrictions. some examples are ets family tfbss, which we detected concentrated around the region immediately upstream of the tss in human and mouse clusters, and nrf- binding sites, which we found in human and mouse clusters, here too in all cases in the region immediately upstream of the tss. importantly, not all motifs were enriched in the upstream regions: zfx motifs, for example, were found to be enriched in regions roughly between to bases downstream of the tss in clusters. standard motif over-representation analyses are likely to miss such regulatory patterns, as this region is usually not included in the input. region with local enrichment of ets domain tf binding motifs were found in different sets of sequences. the species, and the tissues and cell types associated with the promoters in which the motif was found are listed, as well as the regions of local enrichment. tfbss in locally enriched regions tend to have high evolutionary conservation another illustration of the validity of our findings is the relatively high evolutionary conservation of tfbss present within enriched regions, compared to nucleotides within predicted tfbss at the same position relative to tsss of randomly sampled genes. among the enriched regions containing ≥  tfbss that were detected in the human gnf geneatlas sets, contain tfbss that are on average more conserved than randomly picked positions, while only contain sites that are on average less conserved. similarly, for mouse datasets, we found out of regions with highly conserved sites, and only with low conservation. in both human and mouse, the vast majority of the peaks with high conservation were located proximal to the tss. among the relatively more conserved sites are the predicted sites for yy in genes specifically expressed in mouse thymus. the bases corresponding to the yy sites predicted in the region - to + have an average phastcons score of. evolutionary conservation was not only found on the level of predicted tfbss, but also on the level of predicted regions of local enrichment, as illustrated in figure  for promoters of testis-specific genes. for the sites of pwms representing rfx rfx tfs in general, a creb motif, and myb family tfs, similar regions of enrichment were predicted in human and mouse sequences. locally enriched regions tend to be enriched for weak tfbss recently a number of studies, including computational studies using thermodynamic models, have illustrated the importance of weak tfbss in the regulation of transcription. although pwm scores are believed to be correlated with the tf-dna binding energy, the sheer number of sites having only moderately high pwm scores prohibits most computational methods from analyzing them in more detail. one advantage our approach offers is that we can restrict the region of interest to the region preferred by the stronger tfbss, and evaluate if weaker binding sites show a similar local enrichment. our analysis of weak tfbss in the enriched regions found in the human clusters shows that out of regions show a significant enrichment of weaker tfbss. similarly, in mouse datasets, out of regions were enriched for weak tfbss. this result does not include the presence of stronger tfbss, nor weaker sites overlapping with them. for example, promoters of genes expressed in human liver are enriched for hnf tfbss in the region - to + our analysis indicates that in addition to the strong hnf sites, there are weak tfbss predicted in the region of local enrichment in the promoters in this set of genes,. application on sets of genes with differential expression upon tlr stimulation the regulatory pathways controlling toll-like receptor signalling have been relatively well studied, and main regulators such as nf-κb and irfs have been identified. however, it is likely that promoters of genes that are induced or repressed after tlr stimulation contain additional regulatory motifs that are still to be identified. from gene expression data taken from mouse dcs stimulated with tlr ligands, we prepared clusters of co-expressed genes. applying our method to these clusters, we could detect locally enriched regions for the sites of creb-binding tfs, nf-κb, and irfs in clusters of transiently induced genes with peak expression around and 6 hrs after stimulation, respectively. these motifs could also be found by standard motif over-representation analysis. however, our approach also detected a local enrichment for stat binding sites in a cluster with induced expression peaking around 6 hours after stimulation, which could not be detected by the standard approach. in addition, we could make a number of findings that might offer hypotheses for further wet-lab experiments. one is the enrichment of tata box motifs only in the set of promoters corresponding to genes with the fastest transient activation, which suggests that the tata box plays a role in the rapid induction of these genes. finally, a number of enriched regions downstream of the tss were predicted. these include regions enriched in sites for hmx, rfx sites, mtf, zfp and e-box motifs.  although various studies have illustrated that in eukaryotic genomes cis-regulatory motifs can be positioned several kilobases or even megabases away from their target genes, it has also been reported that in a number of cases tfbss show a tendency to be present at a more or less fixed distance with regard to the tss. nevertheless, in general, no clear positional preferences have been described for most regulatory motifs, even though such information could be extremely useful for their prediction. here, we present a novel method for detecting locally enriched tfbss in the regulatory regions of sets of co-regulated genes. our approach is based on a non-parametric approach for sample density estimation, with adjustments which allow it to detect sequence regions that have a significant local enrichment in tfbss, on a 1 bp resolution. an online tool of our approach, which we call locamo finder, is available at http://sysimm.ifrec.osaka-u.ac.jp/tfbs/locamo/. our approach evaluates both positioning and enrichment of tfbss simultaneously, using a set of control sequences as a reference. we implemented our method for the estimation of significance of enrichment in a way that takes into account gc content profiles of the input set of promoter sequences. evaluation of significance is done against randomly selected promoters that have gc content profiles similar to the input sequences. we showed how our approach has a high accuracy compared to other methods and measures for local and global motif enrichment. we applied our method on a large number of sets of genes with tissue- or cell type-specific expression, as well as on a number of sets of genes with similar expression profiles after tlr ligand stimulation in mouse dcs. detected regions of local enrichment of tfbss are supported by known regulatory interactions reported in literature, as well as by the observation that several regulatory motifs are found to be locally enriched in similar regions in different sets of promoters. in addition, detected regions tend to contain tfbss with higher evolutionary conservation than expected, and they also tend to be enriched for weak tfbss. together, these results illustrated the usefulness and validity of our approach. approaches for finding local enrichment of tfbss heavily rely on tss annotations. recent studies have shown that a large fraction of genes have several tsss, and that different tsss might allow different amounts of variation in the bases from which transcription is initiated. future approaches that incorporate such features are likely to present more biological insights into the relationship between tfbs positioning and transcription initiation and tss variability. abbreviations dc: dendritic cell; ori: over-representation index; pwm: position weight matrix; rmsd: root-mean-square deviation; sd: standard deviation; tf: transcription factor; tfbs: transcription factor binding site; tlr: toll-like receptor; tss: transcription start site; zoops: zero or one occurrence per sequence competing interests the authors declare that they have no competing interests. authors’ contributions av conceived of the study and performed the bioinformatics analysis, and prepared the manuscript. yk, st, and kma contributed to statistical analysis, discussion and interpretation of results and with drafting the manuscript. sa and dms supervised the project and helped with discussion and interpretation of results and with drafting the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 supplementary text and figures. a file containing supplementary material, including detailed description of methods and results, as well as supplementary figures and tables. click here for file additional file 2 supplementary tables s s and s a spreadsheet file containing predicted regions of local enrichment and globally enriched tfbs motifs for human and mouse gnf geneatlas datasets, and for the sets of promoters obtained from tlr-stimulated mouse dc cells. click here for file acknowledgements the authors would like to thank all members of the laboratory of systems immunology and the laboratory of host defense for helpful discussions and advice; dr k nakai for helpful suggestions; m ogawa, e kamada and m kageyama for secretarial assistance. we also thank the reviewers for their suggestions, which improved the paper considerably. funding this work was supported by the special coordination funds of the japanese ministry of education, culture, sports, science and technology, and the ministry of health, labour and welfare in japan, the japan society for the promotion of science through funding program for world-leading innovative r&d on science and technology, and by a kakenhi grant-in-aid for scientific research from the japan society for the promotion of science.
the method makes use of a parzen window-based approach for scoring local enrichment, and during evaluation of significance it takes into account gc content of sequences. an online application of our approach, called locamo finder, is available at http://sysimm.ifrec.osaka-u.ac.jp/tfbs/locamo/. we have developed a novel method for the detection of regulatory motifs based on their local over-representation in sets of regulatory regions. one is that binding sites for regulatory proteins are often not significantly over-represented in the set of dna sequences of interest, because of high levels of false positive predictions, and because of positional restrictions on functional binding sites with regard to the transcription start site. we show that the accuracy of our method compares favourably to that of other methods, and that our method is capable of detecting not only generally over-represented regulatory motifs, but also locally over-represented motifs that are often missed by standard motif detection approaches. our approach can be used to suggest testable hypotheses for wet-lab experiments. identification of cis- and trans-acting factors regulating gene expression remains an important problem in biology. bioinformatics analyses of regulatory regions are hampered by several difficulties.
 identification of cis- and trans-acting factors regulating gene expression remains an important problem in biology. bioinformatics analyses of regulatory regions are hampered by several difficulties. one is that binding sites for regulatory proteins are often not significantly over-represented in the set of dna sequences of interest, because of high levels of false positive predictions, and because of positional restrictions on functional binding sites with regard to the transcription start site.  we have developed a novel method for the detection of regulatory motifs based on their local over-representation in sets of regulatory regions. the method makes use of a parzen window-based approach for scoring local enrichment, and during evaluation of significance it takes into account gc content of sequences. we show that the accuracy of our method compares favourably to that of other methods, and that our method is capable of detecting not only generally over-represented regulatory motifs, but also locally over-represented motifs that are often missed by standard motif detection approaches. using a number of examples we illustrate the validity of our approach and suggest applications, such as the analysis of weaker binding sites.  our approach can be used to suggest testable hypotheses for wet-lab experiments. it has potential for future analyses, such as the prediction of weaker binding sites. an online application of our approach, called locamo finder, is available at http://sysimm.ifrec.osaka-u.ac.jp/tfbs/locamo/. regulation of transcriptionpromoter sequencetranscription factor binding siteparzen window 