BACKGROUND
transcription factors  are important components of gene transcriptional regulation.
the binding of a tf to a specific locus can be developmental-stage or cell-type specific; alternatively, as growing evidence suggests, sometimes a protein binds to a specific locus constitutively, i.e., in all the cell types studied so far.
a good example is the ccctc-binding factor .
studies using chromatin immunoprecipitation followed by microarray  or sequencing  showed that, unlike many other tfs/proteins, a portion of the ctcf sites are constitutively bound  <cit>  as illustrated in figure  <dig>  we believe that these constitutive binding sites are likely to have unique or fundamental biological roles.
recently, we carried out a comprehensive analysis of the  <dig> ctcf chip-seq datasets from  <dig> cell lines from the encode  consortium and identified ~ <dig>  ctcf binding sites that were bound in more than 90% of the  <dig> cell lines  <cit>  .
because these constitutive ctcf binding sites were enriched among ctcf-mediated long-range chromatin interactions in k <dig> and mcf <dig> cell lines, we hypothesized that these constitutive ctcf binding sites might play a role in maintaining and/or establishing chromatin structures common to all cell types  <cit> .
thus, we see value in locating constitutive binding sites for other dna binding proteins as a possible window into highly conserved biological processes.
many tfs bind to dna directly and have well-defined motif models.
for such tfs, binding sites may be located by scanning their chip-seq peak sequences with motif models like position weight matrices   <cit> .
using chip-seq data for the same tf from a number of cell lines, one would consider a binding site constitutive if it were found in a sufficiently high proportion of the available cell lines/types.
we refer to this as the motif-based approach.
not every transcription factor has a known binding motif, however.
among the ~ <dig>  sequence-specific dna-binding transcription factors in the human genome  <cit> , only about 10-20% of them have known binding motifs  <cit> .
thus, while the motif-based approach should work well for factors with well-defined pwms such as ctcf, it will fail for tfs lacking reliable pwms or for proteins that do not always bind to their target dna sites directly.
a simple alternative approach, one that accommodates tfs that binds indirectly or that lack a well-defined pwm, is to divide the genome into fixed width bins and count the number of peak centers from chip-seq that fall into each bin, e.g.,  <cit>  and  <cit> .
bins containing peak centers for a sufficiently high proportion of the available cell lines/types are declared as constitutive.
although this binning method is simple, intuitive and commonly used in genome analysis, it suffers from several drawbacks, including a boundary effect where which of two adjacent bins contains a peak center may be ambiguous.
our t-kde approach is based on the following idea.
if a particular genomic locus is bound by a protein in all available cell lines, then the centers of all chip-seq peaks, one from each cell line, at that locus should be within close proximity .
we aimed to systematically identify such sites from chip-seq experiments that target the same protein in multiple cell lines by simultaneously analyzing peak centers across all the experiments.
our goal is distinct from peak calling in a single chip-seq experiment.
a chip-seq peak is a genomic region  enriched with sequence reads and identified using a peak-calling algorithm, e.g.
<cit> .
various peak calling algorithms find genomic regions enriched for a binding signal in a variety of data types.
these include hidden markov model -based peak calling algorithms for chip-seq data  <cit> , for chip-chip data  <cit> , and for medip-seq data  <cit> .
all identify peaks by modeling emission and transition probabilities using multiple states and exploiting distinct signal signatures in different states.
despite the distinct goals, another approach for detecting constitutive protein binding sites might be to apply existing peak-calling tools to the original chip-seq reads from multiple cell lines simultaneously, expecting that constitutive binding sites will exhibit especially high peaks.
such an approach has several drawbacks.
first, bam files from individual chip-seq experiments can be very large, so that combining and processing bam files from tens or hundreds of experiments together will be computationally intensive.
secondly, combining read counts from multiple data sets where some binding occurs at loci common across many data sets and some binding occurs at loci specific to particular data sets will introduce unusual patterns variation in reads counts that could bias estimation of background rates.
for tools that require estimation of background models, this feature may compromise their ability to reliably detect constitutive binding sites.
finally, the definition of constitutive binding site in terms of binding in most cell lines to the same site does not directly translate to a criterion based on peak height in a combined bam file.
consequently, declaring a constitutive peak seems to require mapping all reads under each detected peak back to their original bam files – an additional computational burden.
in this paper, we propose an effective and efficient alternative to binning for locating binding sites for tfs that may bind directly or indirectly.
like binning, it uses peak centers from chip-seq as input data.
our algorithm, t-kde, identifies binding site locations by combining a kernel density estimator  with a binary range tree.
kernel density estimation, also known as the parzen window method, is an unsupervised and non-parametric technique for estimating a continuous probability density function from sample data  <cit> .
because kdes can converge asymptotically to any density function  <cit> , they are widely used and have been applied to many genomic problems such as chip-seq peak calling  <cit> , analyzing nucleosome positioning  <cit>  and detecting transcription factor binding motifs based on their over-representation in regulatory regions  <cit> .
in this paper, we use a kde to find those genomic regions that contain the highest density of chip-seq peak centers from multiple cell lines/types for a given tf.
use of a binary range tree in conjunction with kernel density estimation enhances t-kde’s speed.
a binary range tree is a helpful algorithm for many applications involving range or nearest neighbor searches, indexing and clustering  <cit> ; we use it to recursively subdivide the set of peak centers into subgroups that allow efficient density estimation and mode finding.
using information on the location of peak centers from  <dig> ctcf chip-seq datasets from the encode project, we compared t-kde to both the motif-based approach and the binning approach.
t-kde outperformed the binning approach and was competitive with the motif-based approach.
more than 90% of the t-kde-declared constitutive ctcf binding sites were within  <dig> base pairs  from the nearest motif-declared constitutive ctcf binding sites  — indicating that t-kde is highly accurate.
in addition, t-kde also identified additional constitutive ctcf binding sites that the motif-based approach failed to find due to lack of apparent motif sites in the chip-seq peaks.
we also applied t-kde to  <dig> other proteins for which replicate chip-seq datasets were available in six or more cell lines and found that the number of constitutive binding sites varied from less than a hundred to tens of thousands.
gene ontology  analysis of the genes with constitutive binding sites in their promoters suggests that constitutive binding sites for several of the proteins are biologically meaningful.
methods
data
we downloaded data on chip-seq peaks for  <dig> transcription factors  from the encode portal at the ucsc genome browser  <cit> .
for each chip-seq peak, we calculated the location of the peak center as half the sum of the start and end coordinates for the peak, and we used these locations for subsequent analysis.
location of constitutive ctcf binding sites via motif model: our “gold standard”
for each of  <dig> ctcf chip-seq datasets with at least one replicate, we extended/trimmed each peak to  <dig> bp in length from its center.
we then used a custom python code to extract the sequences from the grch <dig> assembly stored locally.
next, we predicted the locations of the ctcf binding sites in the sequences using the gadem software  <cit>  with a ctcf position weight matrix  derived previously  <cit>  .
we declared a subsequence a ctcf binding site when its pwm score exceeded the score corresponding to the p-value cutoff of  <dig> .
when more than one ctcf site was found in the sequence for a single peak, only the highest scoring site  was retained for that peak.
when a ctcf binding site was found in two or more replicate datasets representing a single cell line, the site was declared present in that cell line.
a ctcf binding site was considered a constitutive binding site when the same motif site was present in more than 90% of the  <dig> cell lines.
we used the center of the motif site as the location of the motif-based binding site.
identification of constitutive binding sites via binning
we divided each chromosome of the human genome into bins of equal size beginning at the centromere and proceeding outward along each arm .
the center of any bin containing peak centers from at least  <dig> replicate datasets from the same cell line was declared a binding site location as before, and those bins containing peak centers for more than 90% of the cell lines were declared constitutive.
we examined this binning procedure with various bin sizes ranging from  <dig> to  <dig> bp.
identification of constitutive binding sites via t-kde
overview
for chip-seq datasets from multiple cell lines, t-kde identifies genomic regions where peak centers occur.
t-kde starts by partitioning the locations of peak centers into subsets, called terminal nodes, using a binary range tree algorithm .
for each terminal node, t-kde uses kernel density estimation to estimate a probability density using all peak centers in the node, the relative frequency with which a peak center will occur near each location along the portion of the genome spanned by the terminal node .
in step  <dig>  for each node, t-kde finds all local maxima and minima of the estimated probability density function and uses them to define modal regions.
the location of each local maximum is taken as a “binding site” location.
in this analysis, we required that a modal region contained peak centers from at least two data sets from the same cell line.
this requirement is meant to reduce false positives .
a binding site is declared constitutive when its modal region contains peak centers from at least two replicate datasets per cell line for more than 90%  of the available cell lines.
below, we describe each step in more details; the algorithm is provided in additional file  <dig>
binary range tree
a binary range tree is an algorithm that produces a structure with all data points stored in the leaves  of the tree for efficient data retrieval and manipulation  <cit> .
in our application, we construct a separate range tree for each chromosome.
initially, all peak centers on the chromosome  are ordered from the smallest to the largest according their genomic locations and placed in the top node.
then, the midrange  is used to partition the peak centers into two sub-nodes: the left sub-node contains peak centers whose locations are less than the mid range whereas the right sub-node contains peak centers whose locations equal or exceed the midrange.
this process continues recursively within each sub-node until a stopping criterion is satisfied.
in our case, a sub-node becomes a terminal node when further partitioning it would result in one or two of its children nodes containing peak centers for fewer than 90% of available cell lines.
although each terminal node in our tree contains peak centers from at least 90% of the cell lines, each terminal node may contain zero, one, or more constitutive binding sites as determined by the subsequent kde analysis and mode finding.
kernel density estimation
kernel density estimation provides a way to estimate the probability density function of a random variable without assuming a particular parametric form  <cit> .
for n independent samples {x <dig>  x <dig>  …, x
n
} drawn from the same unknown distribution with density f, a kernel density estimate of f at any point x, f^x, is given by:
f^x=1n∑i=1n1nkx-xih,
where h represents the bandwidth, a user-defined tuning parameter that controls the smoothness of the resulting estimate.
the kernel k is a symmetric  function that integrates to one, i.e., ∫ kdx =  <dig>  the kernel function serves to smear the probability mass of each data point across a local region.
t-kde uses the gaussian kernel, the density function for a gaussian random variable with mean zero and variance one, defined as:
kx=2π-1/2exp-x2/ <dig>
with this kernel, each term in the sum of equation  is a gaussian density with mean xi and standard deviation h. thus, equation  states that the estimate f^x at any location x is formed by averaging contributions from gaussian densities with standard deviation h and means at the observed peak centers.
the basic operations of kernel density estimation used by t-kde have been modified directly from the kde toolbox for matlab  <cit> .
mode finding in gaussian mixture models
to find local maxima and minima of the estimated density function, we adapted a fixed-point iterative search scheme  <cit> .
our kernel density estimate is an equally weighted mixture of gaussian densities where the mean of each component is an observed peak center.
such a gaussian mixture has, at most, as many local maxima as it has components.
if peak centers are far apart relative to the bandwidth, each peak center will yield a local maximum.
if peak centers are close relative to bandwidth, a local maximum must be between their smallest and largest smallest locations.
thus, within each terminal node, we can use a “hill-climbing” algorithm starting from every peak center to locate all the local maxima and minima.
once we find a location whose gradient is zero using newton’s method, we use a second derivative test to determine whether it is a maximum or a minimum.
modal regions are defined as extending from the observed peak center farthest to the left of the local maximum but no farther than the next local minimum to the similarly delimited observed peak center farthest to the right.
gene ontology
we used david  <cit>  to analyze gene ontology .
we assigned a constitutive binding site to a gene if the site was located within ±5kb from the gene’s transcription start site using the ucsc refgenes model .
all unique genes that were within the distance were included in the go analysis.
RESULTS
utility of the binary range tree
without initial data partition using the binary range tree, kde analysis and mode finding on even a single chromosome is computationally prohibitive; estimating the density, rather than finding the local maxima/minima, is the bottleneck.
for the ctcf datasets, analysis of chromosome  <dig> took less than half an hour with the binary range tree compared to more than  <dig> days without it .
the locations of sites declared constitutive using kde with and without the binary range tree were nearly identical .
bandwidth and bin width selection
the bandwidth  of a kernel in kde estimation and the bin width in the binning approach are tuning parameters whose choice influences each method’s performance in locating binding sites .
thus, selection of an appropriate bandwidth or bin width is crucial to accurate identification of binding sites.
we systematically tested several different choices for their performance in identifying binding sites using the  <dig> ctcf chip-seq datasets .
we used both binning and t-kde to identify binding site locations chromosome by chromosome, declaring a binding site constitutive if the corresponding bin or modal region contained peak centers from at least two replicate datasets per cell line for more than 90% of available cell lines.
both the total number of binding sites and the number of constitutive binding sites identified depended on the method  and on the value of the corresponding tuning parameter .
for t-kde, as expected, the total number of modal regions decreased as the bandwidth increased since a larger bandwidth results fewer but broader modal regions; however, the number of constitutive binding sites remains relatively unchanged.
for binning, the number of total binding sites also decreased with increasing bin width because the number of bins, hence, the number of possible binding sites, decreased with increasing bin width, however, the number of constitutive binding sites increased with bin width because wider bins accumulate adjacent binding sites into the same bin, wrongly declaring several non-constitutive sites amalgamated into the same bin as constitutive.
applying the motif-based approach to the same  <dig> ctcf chip-seq data sets with the same criteria  identified  <dig>  constitutive ctcf binding sites .
we regarded those motif-based constitutive ctcf biding sites as an “alloyed gold standard”.
we have high confidence in a ctcf binding site identified by the motif-based approach because binding at the exact same motif location is detected in more than 90% of cell lines.
on the other hand, the motif-based approach is imperfect as it may fail to identify low affinity or indirect binding sites.
the motif-based approach could also overlook constitutive sites if the length of peak sequence scanned  is too short to cover the actual binding site.
to compare the locations of constitutive ctcf binding sites from t-kde and from binning to the locations of our motif-based constitutive ctcf binding sites.
we plotted the proportion of constitutive binding sites identified by each method that are less than distance d from the nearest motif-based constitutive ctcf binding site as a function of distance d : binning in figure 3).
for t-kde with bandwidths smaller than  <dig> bp, all ctcf binding sites declared constitutive are within  <dig> bp of their nearest motif-based constitutive ctcf binding sites.
for a bandwidth of  <dig> bp, more than 90% of the t-kde-declared constitutive ctcf binding sites are within  <dig> bp of the nearest motif-based constitutive ctcf binding sites and nearly all are within  <dig> bp.
for bandwidths exceeding  <dig> bp, performance deteriorates though roughly 90% of the t-kde-declared constitutive binding sites are still within  <dig> bp from their nearest motif-based counterpart.
the results from table  <dig> and figure  <dig> strongly suggest that changing the bandwidth with t-kde has little impact on the number of constitutive binding sites identified but a greater impact on their locations.
on the other hand, changing the bin width with the binning approach has an impact on both the number of constitutive binding sites identified and on their locations.
our results also suggest that, for ctcf, a bandwidth near  <dig> bp and a bin width near  <dig> bp may be the optimal values for t-kde and for the binning method, respectively.
although derived from ctcf comparisons, we believe these choices of bandwidth or bin width should be applicable to other factors whose chip-seq peak length distributions are similar to those of ctcf.
comparing figure  <dig> and  <dig> also reveals that the accuracy of t-kde for locating constitutive binding sites is generally far superior to that of the binning approach.
in particular, the optimal bandwidth of  <dig> bp was more accurate in locating constitutive binding sites than the optimal bin width of  <dig> bp.
consequently, for our remaining analyses, we focus on t-kde using a bandwidth of  <dig> bp.
t-kde versus binning
using our motif-based constitutive ctcf binding sites as the reference, we believe that the optimal bandwidth for t-kde is  <dig> bp and the optimal bin size for binning method is  <dig> bp.
we plotted the proportion of constitutive binding sites identified by each method that are less than distance from the nearest motif-based constitutive ctcf binding site as a function of distance d .
as shown in figure  <dig>  t-kde is much more accurate in locating constitutive ctcf binding sites than binning method.
constitutive sites found by t-kde but not by the motif-based approach
only  <dig> of the  <dig>  motif-based constitutive ctcf binding sites were farther than  <dig> bp from the nearest constitutive ctcf binding sites identified by t-kde.
furthermore, t-kde declared an additional  <dig>  ctcf binding sites as constitutive that the motif-based approach missed.
among those  <dig>  sites, the motif-based approach failed to detect  <dig> because no sub-sequence in the corresponding peaks reached the motif significance cutoff.
.
the remaining additional constitutive sites found by t-kde were found by the motif-based approach in a majority of cell lines but not in enough cell lines to reach the required 90%.
when the true binding sites are not located near the center of some peaks and/or the peak sequences used in motif scan are not long enough to cover the actual motif, a motif-based approach would miss the site.
t-dke, however, is unaffected by these issues.
because it uses peak centers from all cell lines to identify the center of mass of each modal region as the binding site, some misalignment or displacement among chip-seq peaks is tolerated.
thus, t-kde is capable of identifying constitutive binding sites that are bound by a protein either directly or indirectly.
motif discovery using gadem  <cit>  on the  <dig> constitutive sites identified by t-kde where no canonical ctcf motif was found yielded two new motif variants  of the canonical core ctcf motif.
we named these motifs as core variant motifs  <dig> and  <dig> , respectively.
cv <dig> was the dominant motif found in ~65% of the  <dig> sequences whereas cv <dig> was found in ~35% of the sequences.
compared to the canonical core ctcf motif, these motif variants lacked information at either the 5′- or the 3′-end.
this feature may explain their lack of motif significance when compared to the canonical core ctcf motif.
both variants are likely authentic, as they are highly centrally distributed along the  <dig> bp peak sequences .
it is also likely that such variant motifs are not limited to those  <dig> constitutive ctcf sites.
the two new motifs along with those discovered recently  <cit>  add to the complexity of the ctcf code.
analysis of constitutive binding sites for  <dig> factors
we applied t-kde using a bandwidth of  <dig> bp to  <dig> factors with chip-seq data available from encode for multiple cell lines with replicates.
the number of declared binding sites  among these tfs ranged from  <dig>  to over  <dig> , and the number constitutive binding sites ranged from a few to over  <dig>  .
constitutive binding sites identified for tfs that were studied in fewer than  <dig> cell lines, especially when the number of constitutive binding sites identified is relatively small, have a high likelihood of being false positives.
as data from additional cell lines becomes available, some sites now declared constitutive could fail to meet the necessary criterion.
thus, we focus attention on the possible biological roles of constitutive binding sites for tfs with more than  <dig>  declared constitutive binding sites based on  <dig> or more cell lines.
besides ctcf and rad <dig>  pol ii exhibited a large number of constitutive binding sites , although they represented only a small proportion of the total pol ii binding sites.
gene ontology  analysis of the genes containing constitutive pol ii binding sites in their proximal promoters suggests that those pol ii target genes are highly enriched with biological processes such as metabolism, biosynthesis and cell cycle .
similarly, genes with constitutive binding sites for other tfs are highly enriched in certain biological processes.
these include gabp, nrsf, taf <dig>  etc.
.
together, those results suggest that binding sites declared constitutive by t-kde are connected to important biological processes.
p
discussion
binding sites that are occupied by a protein regardless of the cell or tissue type seem likely to have a distinct role compared to binding sites for the same protein that are occupied more selectively – the constitutive nature of the binding should signify something of fundamental import.
our earlier work using motif-based analysis found that constitutive ctcf binding sites, especially those near rad <dig> sites, are highly enriched in ctcf-mediated chromatin interactions  <cit>  and those interactions are predominately within topological domains, not between them  <cit> .
consequently, we hypothesized that the constitutive ctcf binding sites may be involved in maintaining and/or establishing chromatin structures that are common among most human cell types  <cit> .
those earlier findings indicate to us that constitutive binding sites for other tfs may have unique biological roles.
the encode consortium has generated more than  <dig>  chip-seq protein-binding datasets for more than  <dig> proteins in multiple cell lines, and the data continue to expand.
discovering the locations and functions the genomic loci that are constitutively bound by each of the proteins is potentially important.
however, computational methods for locating constitutive binding sites when the protein does not bind directly to dna are still lacking.
one challenge is that the chip-seq peak data are low-resolution, and the technology is unable pinpoint exact genomic binding locations.
to fill this gap, we developed an efficient and effective approach, t-kde which takes as input locations of peak centers from multiple chip-seq data sets and returns estimates of the locations of binding sites and declares them constitutive or not.
t-kde combines a binary range tree algorithm, a kernel density estimator, and a mode finding algorithm.
using data on ctcf binding, we found that t-kde was superior at locating constitutive binding sites compared to a naïve approach based on binning and that t-kde performed well compared to the motif-based approach.
for example, all motif-based constitutive ctcf binding sites were included in the constitutive ctcf binding sites identified by t-kde.
furthermore, t-kde identified additional  <dig>  constitutive ctcf binding sites that the motif-based approach failed to detect.
this result highlights a major advantage of t-kde compared to both the motif-based and binning approaches: regardless of whether binding is direct or indirect and whether an adequate motif model is known, t-kde accurately estimates the locations of constitutive binding sites by identifying genomic regions where the centers of chip-seq peaks from multiple datasets lie in close proximity.
accurate binding locations are necessary for subsequent functional analysis and discovery.
we applied t-kde to locate constitutive binding sites, if present, for  <dig> tfs that had replicate chip-seq data sets for at least  <dig> cell lines available from encode, and we used gene ontology analysis to establish possible biological functions for some of those tfs.
kde-based methods different from ours have been applied to chip-seq reads for peak calling  <cit>  and nucleosome positioning  <cit> .
additionally, kde-based method has been applied to motif locations for detection of regions locally enriched with transcription factor binding sites  <cit> .
our goal is different: we use the locations of chip-seq peak centers from multiple cell lines  to infer the location of constitutive binding sites.
in addition, our method has unique features.
our method first recursively partitions the locations of peak centers into subgroups  using a binary range tree algorithm.
the partitioning stops whenever either of the two would-be child nodes contains peak centers from fewer than 90%  of available cell lines.
the kde analysis and subsequent mode finding is carried out on each terminal node, one at a time.
the partitioning guarantees that more than 90% of cell lines are represented in every terminal node; however, a terminal node may still contain zero, one or more constitutive binding sites depending on the spread of the peak centers present — making kde and subsequent mode-finding necessary for localizing modal regions.
binding site locations are declared at local maxima within modal regions.
our use of the binary range tree before applying kde and mode-finding makes our algorithm novel and efficient.
one reviewer suggested an alternative procedure  using the peak-finding algorithm macs  <cit> .
the procedure involves applying macs in its default parameters to a combined bam file from the original chip-seq reads data  from the multiple cell lines.
the peaks with low variation in log  within ± <dig> bp from the macs summit are considered constitutive.
we compared this procedure with t-kde and a binding-based method and showed that t-kde was far superior to this alternative procedure .
although t-kde can be applied to chip-seq data from any number of cell lines, caution must be excised when interpreting a result from only a few cell lines.
because the property of being constitutive requires binding to the same locus in a variety of cell types, the number and diversity  of cell lines/types providing data to the algorithm would be expected to have a strong influence on the biological trustworthiness of any result.
for n peak centers, kde followed by mode-finding has a computational complexity of o  <cit> .
when n is large as in our ctcf dataset , the process becomes computationally prohibitive.
after initial data partitioning by a binary range tree into a set of terminal nodes indexed by i, each with n peak centers, complexity is greatly reduced to   ∑ io.
consequently, t-kde reduces the computational time for ctcf on chromosome  <dig> from days to within an hour.
we envision that parallelization of our t-kde algorithm at the node level would further reduce the computational time.
a potential cost is that partitioning all peak centers onto terminal nodes before the kde analysis and mode finding might destroy a constitutive binding site by splitting it between two adjacent nodes.
this problem appears to arise rarely or not at all as we observed that the performance of t-kde was nearly identical to that of kde omitting the initial partitioning.
we attribute this similarity, in part, to our stopping criterion for partitioning.
generally, the choice of the bandwidth for kde can exhibit a strong influence on the shape of the estimated density: small bandwidths yielding spiky estimates and large bandwidths yielding overly flattened ones.
yet, in our comparisons when locating constitutive ctcf binding sites, bandwidths from  <dig> to  <dig> bp uncovered similar numbers of constitutive ctcf binding sites and the distribution of the distances from t-kde-declared sites to the nearest motif-declared sites did not change much with bandwidth.
we believe that a bandwidth of  <dig> to  <dig> bp may be optimal for most tf binding sites with narrow peaks .
automatic selection of the optimal bandwidth would be desirable, but optimal bandwidth selection based statistical criteria such as the mean integrated squared error  <cit>  did not work well with the ctcf data.
that process, which involved maximizing a “pseudo-likelihood” combined with a leave-one-out cross-validation approach  <cit>  was computationally expensive and selected a large bandwidth of  <dig>  bp that did not locate constitutive binding sites as well as our preferred  <dig> bp bandwidth did.
although designed for identifying constitutive binding sites for a protein using chip-seq data from multiple cell lines, our method could also be used to identify genomic loci that have concentrations of different protein binding sites , and conversely “cold spots”, using multiple protein chip-seq data for the cell line.
CONCLUSIONS
in conclusion, we developed efficient and accurate method, t-kde, to locate constitutive protein binding sites using chip-seq peak centers from multiple cell lines.
t-kde combines a binary range tree algorithm, a non-parametric kernel density estimator, and a mode finding algorithm.
we showed that, for ctcf data, our method is relatively robust to the choice of bandwidth and is highly accurate when compared to the identification of constitutive binding sites through motif analysis.
application of t-kde to  <dig> proteins with chip-seq data from multiple cell lines located substantial numbers of constitutive binding sites for some tfs but almost none for others.
for tfs with large numbers of constitutive binding sites, go analysis suggests that these sites are biological meaningful.
as additional tf binding sites chip-seq datasets become available in more cell lines and for more tfs, our method will prove to be essential for identifying their constitutive binding sites.
availability and requirements
project name: t-kde
project homepage: http://www.niehs.nih.gov/research/resources/software/biostatistics/t-kde/index.cfm
operating system: unix
programming language: matlab
other requirements: n/a
license: this work is made available under the gpl v <dig>
any restrictions to use by non-academics: none
abbreviations
chip-chip: chromatin immunoprecipitation followed by microarray; chip-seq: chromatin immunoprecipitation followed by sequencing; ctcf: ccctc binding factor; gabp: ga-binding protein; go: gene ontology; kde: kernel density estimation; macs: model-based analysis of chip-seq; nrsf: neuron-restrictive silencer transcription factor, also known as rest; pol ii: polymerase ii; rad21: homolog ; taf1: tata box binding protein-associated factor 1; tf: transcription factor.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
ll conceived the study, yl, and ll performed the analyses, and yl, ll, and dmu were all involved in study design, planning analyses, and the interpretation of results.
all authors contributed to writing and revising the manuscript, and all approved the final manuscript for publication.
supplementary material
additional file 1: supplementary tables s1-19
cell lines contributing chip-seq data for each of the  <dig> transcription factors.
click here for file
additional file 2
outlines of algorithms.
t-kde.
click here for file
additional file 3: supplementary figure s1
proportion of tkde-declared versus kde-declared constitutive ctcf binding sites whose distance from nearest motif-based constitutive ctcf binding site on  <dig> chromosomesare less than distance d plotted as a function of d for various bandwidths.separate curves for t-kde with bandwidth of  <dig> bps and for the same density estimation algorithm without the binary range tree pre-processing.
click here for file
additional file 4
additional results comparing various methods.
t-kde: a method for genome-wide identification of constitutive protein binding sites from multiple chip-seq data sets.
click here for file
acknowledgements
we thank liang niu and weichun huang for discussion and xuting wang and grace kissling for critical reading of the manuscript.
we thank the computational biology facility at niehs for computing time and support.
this research was supported by intramural research program of the nih, national institute of environmental health sciences .

 transcription factors are important components of gene transcriptional regulation. the binding of a tf to a specific locus can be developmental-stage or cell-type specific; alternatively, as growing evidence suggests, sometimes a protein binds to a specific locus constitutively, i.e., in all the cell types studied so far. a good example is the ccctc-binding factor. studies using chromatin immunoprecipitation followed by microarray or sequencing showed that, unlike many other tfs/proteins, a portion of the ctcf sites are constitutively bound as illustrated in figure  we believe that these constitutive binding sites are likely to have unique or fundamental biological roles. recently, we carried out a comprehensive analysis of the ctcf chip-seq datasets from cell lines from the encode consortium and identified ~ ctcf binding sites that were bound in more than 90% of the cell lines. because these constitutive ctcf binding sites were enriched among ctcf-mediated long-range chromatin interactions in k and mcf cell lines, we hypothesized that these constitutive ctcf binding sites might play a role in maintaining and/or establishing chromatin structures common to all cell types. thus, we see value in locating constitutive binding sites for other dna binding proteins as a possible window into highly conserved biological processes. many tfs bind to dna directly and have well-defined motif models. for such tfs, binding sites may be located by scanning their chip-seq peak sequences with motif models like position weight matrices. using chip-seq data for the same tf from a number of cell lines, one would consider a binding site constitutive if it were found in a sufficiently high proportion of the available cell lines/types. we refer to this as the motif-based approach. not every transcription factor has a known binding motif, however. among the ~ sequence-specific dna-binding transcription factors in the human genome, only about 10-20% of them have known binding motifs. thus, while the motif-based approach should work well for factors with well-defined pwms such as ctcf, it will fail for tfs lacking reliable pwms or for proteins that do not always bind to their target dna sites directly. a simple alternative approach, one that accommodates tfs that binds indirectly or that lack a well-defined pwm, is to divide the genome into fixed width bins and count the number of peak centers from chip-seq that fall into each bin, e.g., and. bins containing peak centers for a sufficiently high proportion of the available cell lines/types are declared as constitutive. although this binning method is simple, intuitive and commonly used in genome analysis, it suffers from several drawbacks, including a boundary effect where which of two adjacent bins contains a peak center may be ambiguous. our t-kde approach is based on the following idea. if a particular genomic locus is bound by a protein in all available cell lines, then the centers of all chip-seq peaks, one from each cell line, at that locus should be within close proximity. we aimed to systematically identify such sites from chip-seq experiments that target the same protein in multiple cell lines by simultaneously analyzing peak centers across all the experiments. our goal is distinct from peak calling in a single chip-seq experiment. a chip-seq peak is a genomic region enriched with sequence reads and identified using a peak-calling algorithm, e.g.. various peak calling algorithms find genomic regions enriched for a binding signal in a variety of data types. these include hidden markov model -based peak calling algorithms for chip-seq data, for chip-chip data, and for medip-seq data. all identify peaks by modeling emission and transition probabilities using multiple states and exploiting distinct signal signatures in different states. despite the distinct goals, another approach for detecting constitutive protein binding sites might be to apply existing peak-calling tools to the original chip-seq reads from multiple cell lines simultaneously, expecting that constitutive binding sites will exhibit especially high peaks. such an approach has several drawbacks. first, bam files from individual chip-seq experiments can be very large, so that combining and processing bam files from tens or hundreds of experiments together will be computationally intensive. secondly, combining read counts from multiple data sets where some binding occurs at loci common across many data sets and some binding occurs at loci specific to particular data sets will introduce unusual patterns variation in reads counts that could bias estimation of background rates. for tools that require estimation of background models, this feature may compromise their ability to reliably detect constitutive binding sites. finally, the definition of constitutive binding site in terms of binding in most cell lines to the same site does not directly translate to a criterion based on peak height in a combined bam file. consequently, declaring a constitutive peak seems to require mapping all reads under each detected peak back to their original bam files – an additional computational burden. in this paper, we propose an effective and efficient alternative to binning for locating binding sites for tfs that may bind directly or indirectly. like binning, it uses peak centers from chip-seq as input data. our algorithm, t-kde, identifies binding site locations by combining a kernel density estimator with a binary range tree. kernel density estimation, also known as the parzen window method, is an unsupervised and non-parametric technique for estimating a continuous probability density function from sample data. because kdes can converge asymptotically to any density function, they are widely used and have been applied to many genomic problems such as chip-seq peak calling, analyzing nucleosome positioning and detecting transcription factor binding motifs based on their over-representation in regulatory regions. in this paper, we use a kde to find those genomic regions that contain the highest density of chip-seq peak centers from multiple cell lines/types for a given tf. use of a binary range tree in conjunction with kernel density estimation enhances t-kde’s speed. a binary range tree is a helpful algorithm for many applications involving range or nearest neighbor searches, indexing and clustering; we use it to recursively subdivide the set of peak centers into subgroups that allow efficient density estimation and mode finding. using information on the location of peak centers from ctcf chip-seq datasets from the encode project, we compared t-kde to both the motif-based approach and the binning approach. t-kde outperformed the binning approach and was competitive with the motif-based approach. more than 90% of the t-kde-declared constitutive ctcf binding sites were within base pairs from the nearest motif-declared constitutive ctcf binding sites — indicating that t-kde is highly accurate. in addition, t-kde also identified additional constitutive ctcf binding sites that the motif-based approach failed to find due to lack of apparent motif sites in the chip-seq peaks. we also applied t-kde to other proteins for which replicate chip-seq datasets were available in six or more cell lines and found that the number of constitutive binding sites varied from less than a hundred to tens of thousands. gene ontology analysis of the genes with constitutive binding sites in their promoters suggests that constitutive binding sites for several of the proteins are biologically meaningful. methods data we downloaded data on chip-seq peaks for transcription factors from the encode portal at the ucsc genome browser. for each chip-seq peak, we calculated the location of the peak center as half the sum of the start and end coordinates for the peak, and we used these locations for subsequent analysis. location of constitutive ctcf binding sites via motif model: our “gold standard” for each of ctcf chip-seq datasets with at least one replicate, we extended/trimmed each peak to bp in length from its center. we then used a custom python code to extract the sequences from the grch assembly stored locally. next, we predicted the locations of the ctcf binding sites in the sequences using the gadem software with a ctcf position weight matrix derived previously. we declared a subsequence a ctcf binding site when its pwm score exceeded the score corresponding to the p-value cutoff of. when more than one ctcf site was found in the sequence for a single peak, only the highest scoring site was retained for that peak. when a ctcf binding site was found in two or more replicate datasets representing a single cell line, the site was declared present in that cell line. a ctcf binding site was considered a constitutive binding site when the same motif site was present in more than 90% of the cell lines. we used the center of the motif site as the location of the motif-based binding site. identification of constitutive binding sites via binning we divided each chromosome of the human genome into bins of equal size beginning at the centromere and proceeding outward along each arm. the center of any bin containing peak centers from at least replicate datasets from the same cell line was declared a binding site location as before, and those bins containing peak centers for more than 90% of the cell lines were declared constitutive. we examined this binning procedure with various bin sizes ranging from to bp. identification of constitutive binding sites via t-kde overview for chip-seq datasets from multiple cell lines, t-kde identifies genomic regions where peak centers occur. t-kde starts by partitioning the locations of peak centers into subsets, called terminal nodes, using a binary range tree algorithm. for each terminal node, t-kde uses kernel density estimation to estimate a probability density using all peak centers in the node, the relative frequency with which a peak center will occur near each location along the portion of the genome spanned by the terminal node. in step for each node, t-kde finds all local maxima and minima of the estimated probability density function and uses them to define modal regions. the location of each local maximum is taken as a “binding site” location. in this analysis, we required that a modal region contained peak centers from at least two data sets from the same cell line. this requirement is meant to reduce false positives. a binding site is declared constitutive when its modal region contains peak centers from at least two replicate datasets per cell line for more than 90% of the available cell lines. below, we describe each step in more details; the algorithm is provided in additional file binary range tree a binary range tree is an algorithm that produces a structure with all data points stored in the leaves of the tree for efficient data retrieval and manipulation. in our application, we construct a separate range tree for each chromosome. initially, all peak centers on the chromosome are ordered from the smallest to the largest according their genomic locations and placed in the top node. then, the midrange is used to partition the peak centers into two sub-nodes: the left sub-node contains peak centers whose locations are less than the mid range whereas the right sub-node contains peak centers whose locations equal or exceed the midrange. this process continues recursively within each sub-node until a stopping criterion is satisfied. in our case, a sub-node becomes a terminal node when further partitioning it would result in one or two of its children nodes containing peak centers for fewer than 90% of available cell lines. although each terminal node in our tree contains peak centers from at least 90% of the cell lines, each terminal node may contain zero, one, or more constitutive binding sites as determined by the subsequent kde analysis and mode finding. kernel density estimation kernel density estimation provides a way to estimate the probability density function of a random variable without assuming a particular parametric form. for n independent samples {x  x  …, x n } drawn from the same unknown distribution with density f, a kernel density estimate of f at any point x, f^x, is given by: f^x=1n∑i=1n1nkx-xih, where h represents the bandwidth, a user-defined tuning parameter that controls the smoothness of the resulting estimate. the kernel k is a symmetric function that integrates to one, i.e., ∫ kdx =  the kernel function serves to smear the probability mass of each data point across a local region. t-kde uses the gaussian kernel, the density function for a gaussian random variable with mean zero and variance one, defined as: kx=2π-1/2exp-x2/ with this kernel, each term in the sum of equation is a gaussian density with mean xi and standard deviation h. thus, equation states that the estimate f^x at any location x is formed by averaging contributions from gaussian densities with standard deviation h and means at the observed peak centers. the basic operations of kernel density estimation used by t-kde have been modified directly from the kde toolbox for matlab. mode finding in gaussian mixture models to find local maxima and minima of the estimated density function, we adapted a fixed-point iterative search scheme. our kernel density estimate is an equally weighted mixture of gaussian densities where the mean of each component is an observed peak center. such a gaussian mixture has, at most, as many local maxima as it has components. if peak centers are far apart relative to the bandwidth, each peak center will yield a local maximum. if peak centers are close relative to bandwidth, a local maximum must be between their smallest and largest smallest locations. thus, within each terminal node, we can use a “hill-climbing” algorithm starting from every peak center to locate all the local maxima and minima. once we find a location whose gradient is zero using newton’s method, we use a second derivative test to determine whether it is a maximum or a minimum. modal regions are defined as extending from the observed peak center farthest to the left of the local maximum but no farther than the next local minimum to the similarly delimited observed peak center farthest to the right. gene ontology we used david to analyze gene ontology. we assigned a constitutive binding site to a gene if the site was located within ±5kb from the gene’s transcription start site using the ucsc refgenes model. all unique genes that were within the distance were included in the go analysis.  utility of the binary range tree without initial data partition using the binary range tree, kde analysis and mode finding on even a single chromosome is computationally prohibitive; estimating the density, rather than finding the local maxima/minima, is the bottleneck. for the ctcf datasets, analysis of chromosome took less than half an hour with the binary range tree compared to more than days without it. the locations of sites declared constitutive using kde with and without the binary range tree were nearly identical. bandwidth and bin width selection the bandwidth of a kernel in kde estimation and the bin width in the binning approach are tuning parameters whose choice influences each method’s performance in locating binding sites. thus, selection of an appropriate bandwidth or bin width is crucial to accurate identification of binding sites. we systematically tested several different choices for their performance in identifying binding sites using the ctcf chip-seq datasets. we used both binning and t-kde to identify binding site locations chromosome by chromosome, declaring a binding site constitutive if the corresponding bin or modal region contained peak centers from at least two replicate datasets per cell line for more than 90% of available cell lines. both the total number of binding sites and the number of constitutive binding sites identified depended on the method and on the value of the corresponding tuning parameter. for t-kde, as expected, the total number of modal regions decreased as the bandwidth increased since a larger bandwidth results fewer but broader modal regions; however, the number of constitutive binding sites remains relatively unchanged. for binning, the number of total binding sites also decreased with increasing bin width because the number of bins, hence, the number of possible binding sites, decreased with increasing bin width, however, the number of constitutive binding sites increased with bin width because wider bins accumulate adjacent binding sites into the same bin, wrongly declaring several non-constitutive sites amalgamated into the same bin as constitutive. applying the motif-based approach to the same ctcf chip-seq data sets with the same criteria identified constitutive ctcf binding sites. we regarded those motif-based constitutive ctcf biding sites as an “alloyed gold standard”. we have high confidence in a ctcf binding site identified by the motif-based approach because binding at the exact same motif location is detected in more than 90% of cell lines. on the other hand, the motif-based approach is imperfect as it may fail to identify low affinity or indirect binding sites. the motif-based approach could also overlook constitutive sites if the length of peak sequence scanned is too short to cover the actual binding site. to compare the locations of constitutive ctcf binding sites from t-kde and from binning to the locations of our motif-based constitutive ctcf binding sites. we plotted the proportion of constitutive binding sites identified by each method that are less than distance d from the nearest motif-based constitutive ctcf binding site as a function of distance d: binning in figure 3). for t-kde with bandwidths smaller than bp, all ctcf binding sites declared constitutive are within bp of their nearest motif-based constitutive ctcf binding sites. for a bandwidth of bp, more than 90% of the t-kde-declared constitutive ctcf binding sites are within bp of the nearest motif-based constitutive ctcf binding sites and nearly all are within bp. for bandwidths exceeding bp, performance deteriorates though roughly 90% of the t-kde-declared constitutive binding sites are still within bp from their nearest motif-based counterpart. the results from table  and figure  strongly suggest that changing the bandwidth with t-kde has little impact on the number of constitutive binding sites identified but a greater impact on their locations. on the other hand, changing the bin width with the binning approach has an impact on both the number of constitutive binding sites identified and on their locations. our results also suggest that, for ctcf, a bandwidth near bp and a bin width near bp may be the optimal values for t-kde and for the binning method, respectively. although derived from ctcf comparisons, we believe these choices of bandwidth or bin width should be applicable to other factors whose chip-seq peak length distributions are similar to those of ctcf. comparing figure  and also reveals that the accuracy of t-kde for locating constitutive binding sites is generally far superior to that of the binning approach. in particular, the optimal bandwidth of bp was more accurate in locating constitutive binding sites than the optimal bin width of bp. consequently, for our remaining analyses, we focus on t-kde using a bandwidth of bp. t-kde versus binning using our motif-based constitutive ctcf binding sites as the reference, we believe that the optimal bandwidth for t-kde is bp and the optimal bin size for binning method is bp. we plotted the proportion of constitutive binding sites identified by each method that are less than distance from the nearest motif-based constitutive ctcf binding site as a function of distance d. as shown in figure  t-kde is much more accurate in locating constitutive ctcf binding sites than binning method. constitutive sites found by t-kde but not by the motif-based approach only of the motif-based constitutive ctcf binding sites were farther than bp from the nearest constitutive ctcf binding sites identified by t-kde. furthermore, t-kde declared an additional ctcf binding sites as constitutive that the motif-based approach missed. among those sites, the motif-based approach failed to detect because no sub-sequence in the corresponding peaks reached the motif significance cutoff.. the remaining additional constitutive sites found by t-kde were found by the motif-based approach in a majority of cell lines but not in enough cell lines to reach the required 90%. when the true binding sites are not located near the center of some peaks and/or the peak sequences used in motif scan are not long enough to cover the actual motif, a motif-based approach would miss the site. t-dke, however, is unaffected by these issues. because it uses peak centers from all cell lines to identify the center of mass of each modal region as the binding site, some misalignment or displacement among chip-seq peaks is tolerated. thus, t-kde is capable of identifying constitutive binding sites that are bound by a protein either directly or indirectly. motif discovery using gadem on the constitutive sites identified by t-kde where no canonical ctcf motif was found yielded two new motif variants of the canonical core ctcf motif. we named these motifs as core variant motifs and, respectively. cv was the dominant motif found in ~65% of the sequences whereas cv was found in ~35% of the sequences. compared to the canonical core ctcf motif, these motif variants lacked information at either the 5′- or the 3′-end. this feature may explain their lack of motif significance when compared to the canonical core ctcf motif. both variants are likely authentic, as they are highly centrally distributed along the bp peak sequences. it is also likely that such variant motifs are not limited to those constitutive ctcf sites. the two new motifs along with those discovered recently add to the complexity of the ctcf code. analysis of constitutive binding sites for factors we applied t-kde using a bandwidth of bp to factors with chip-seq data available from encode for multiple cell lines with replicates. the number of declared binding sites among these tfs ranged from to over, and the number constitutive binding sites ranged from a few to over. constitutive binding sites identified for tfs that were studied in fewer than cell lines, especially when the number of constitutive binding sites identified is relatively small, have a high likelihood of being false positives. as data from additional cell lines becomes available, some sites now declared constitutive could fail to meet the necessary criterion. thus, we focus attention on the possible biological roles of constitutive binding sites for tfs with more than declared constitutive binding sites based on or more cell lines. besides ctcf and rad pol ii exhibited a large number of constitutive binding sites, although they represented only a small proportion of the total pol ii binding sites. gene ontology analysis of the genes containing constitutive pol ii binding sites in their proximal promoters suggests that those pol ii target genes are highly enriched with biological processes such as metabolism, biosynthesis and cell cycle. similarly, genes with constitutive binding sites for other tfs are highly enriched in certain biological processes. these include gabp, nrsf, taf etc.. together, those results suggest that binding sites declared constitutive by t-kde are connected to important biological processes. p discussion binding sites that are occupied by a protein regardless of the cell or tissue type seem likely to have a distinct role compared to binding sites for the same protein that are occupied more selectively – the constitutive nature of the binding should signify something of fundamental import. our earlier work using motif-based analysis found that constitutive ctcf binding sites, especially those near rad sites, are highly enriched in ctcf-mediated chromatin interactions and those interactions are predominately within topological domains, not between them. consequently, we hypothesized that the constitutive ctcf binding sites may be involved in maintaining and/or establishing chromatin structures that are common among most human cell types. those earlier findings indicate to us that constitutive binding sites for other tfs may have unique biological roles. the encode consortium has generated more than chip-seq protein-binding datasets for more than proteins in multiple cell lines, and the data continue to expand. discovering the locations and functions the genomic loci that are constitutively bound by each of the proteins is potentially important. however, computational methods for locating constitutive binding sites when the protein does not bind directly to dna are still lacking. one challenge is that the chip-seq peak data are low-resolution, and the technology is unable pinpoint exact genomic binding locations. to fill this gap, we developed an efficient and effective approach, t-kde which takes as input locations of peak centers from multiple chip-seq data sets and returns estimates of the locations of binding sites and declares them constitutive or not. t-kde combines a binary range tree algorithm, a kernel density estimator, and a mode finding algorithm. using data on ctcf binding, we found that t-kde was superior at locating constitutive binding sites compared to a naïve approach based on binning and that t-kde performed well compared to the motif-based approach. for example, all motif-based constitutive ctcf binding sites were included in the constitutive ctcf binding sites identified by t-kde. furthermore, t-kde identified additional constitutive ctcf binding sites that the motif-based approach failed to detect. this result highlights a major advantage of t-kde compared to both the motif-based and binning approaches: regardless of whether binding is direct or indirect and whether an adequate motif model is known, t-kde accurately estimates the locations of constitutive binding sites by identifying genomic regions where the centers of chip-seq peaks from multiple datasets lie in close proximity. accurate binding locations are necessary for subsequent functional analysis and discovery. we applied t-kde to locate constitutive binding sites, if present, for tfs that had replicate chip-seq data sets for at least cell lines available from encode, and we used gene ontology analysis to establish possible biological functions for some of those tfs. kde-based methods different from ours have been applied to chip-seq reads for peak calling and nucleosome positioning. additionally, kde-based method has been applied to motif locations for detection of regions locally enriched with transcription factor binding sites. our goal is different: we use the locations of chip-seq peak centers from multiple cell lines to infer the location of constitutive binding sites. in addition, our method has unique features. our method first recursively partitions the locations of peak centers into subgroups using a binary range tree algorithm. the partitioning stops whenever either of the two would-be child nodes contains peak centers from fewer than 90% of available cell lines. the kde analysis and subsequent mode finding is carried out on each terminal node, one at a time. the partitioning guarantees that more than 90% of cell lines are represented in every terminal node; however, a terminal node may still contain zero, one or more constitutive binding sites depending on the spread of the peak centers present — making kde and subsequent mode-finding necessary for localizing modal regions. binding site locations are declared at local maxima within modal regions. our use of the binary range tree before applying kde and mode-finding makes our algorithm novel and efficient. one reviewer suggested an alternative procedure using the peak-finding algorithm macs. the procedure involves applying macs in its default parameters to a combined bam file from the original chip-seq reads data from the multiple cell lines. the peaks with low variation in log within ± bp from the macs summit are considered constitutive. we compared this procedure with t-kde and a binding-based method and showed that t-kde was far superior to this alternative procedure. although t-kde can be applied to chip-seq data from any number of cell lines, caution must be excised when interpreting a result from only a few cell lines. because the property of being constitutive requires binding to the same locus in a variety of cell types, the number and diversity of cell lines/types providing data to the algorithm would be expected to have a strong influence on the biological trustworthiness of any result. for n peak centers, kde followed by mode-finding has a computational complexity of o. when n is large as in our ctcf dataset, the process becomes computationally prohibitive. after initial data partitioning by a binary range tree into a set of terminal nodes indexed by i, each with n peak centers, complexity is greatly reduced to   ∑ io. consequently, t-kde reduces the computational time for ctcf on chromosome from days to within an hour. we envision that parallelization of our t-kde algorithm at the node level would further reduce the computational time. a potential cost is that partitioning all peak centers onto terminal nodes before the kde analysis and mode finding might destroy a constitutive binding site by splitting it between two adjacent nodes. this problem appears to arise rarely or not at all as we observed that the performance of t-kde was nearly identical to that of kde omitting the initial partitioning. we attribute this similarity, in part, to our stopping criterion for partitioning. generally, the choice of the bandwidth for kde can exhibit a strong influence on the shape of the estimated density: small bandwidths yielding spiky estimates and large bandwidths yielding overly flattened ones. yet, in our comparisons when locating constitutive ctcf binding sites, bandwidths from to bp uncovered similar numbers of constitutive ctcf binding sites and the distribution of the distances from t-kde-declared sites to the nearest motif-declared sites did not change much with bandwidth. we believe that a bandwidth of to bp may be optimal for most tf binding sites with narrow peaks. automatic selection of the optimal bandwidth would be desirable, but optimal bandwidth selection based statistical criteria such as the mean integrated squared error did not work well with the ctcf data. that process, which involved maximizing a “pseudo-likelihood” combined with a leave-one-out cross-validation approach was computationally expensive and selected a large bandwidth of bp that did not locate constitutive binding sites as well as our preferred bp bandwidth did. although designed for identifying constitutive binding sites for a protein using chip-seq data from multiple cell lines, our method could also be used to identify genomic loci that have concentrations of different protein binding sites, and conversely “cold spots”, using multiple protein chip-seq data for the cell line.  in conclusion, we developed efficient and accurate method, t-kde, to locate constitutive protein binding sites using chip-seq peak centers from multiple cell lines. t-kde combines a binary range tree algorithm, a non-parametric kernel density estimator, and a mode finding algorithm. we showed that, for ctcf data, our method is relatively robust to the choice of bandwidth and is highly accurate when compared to the identification of constitutive binding sites through motif analysis. application of t-kde to proteins with chip-seq data from multiple cell lines located substantial numbers of constitutive binding sites for some tfs but almost none for others. for tfs with large numbers of constitutive binding sites, go analysis suggests that these sites are biological meaningful. as additional tf binding sites chip-seq datasets become available in more cell lines and for more tfs, our method will prove to be essential for identifying their constitutive binding sites. availability and requirements project name: t-kde project homepage: http://www.niehs.nih.gov/research/resources/software/biostatistics/t-kde/index.cfm operating system: unix programming language: matlab other requirements: n/a license: this work is made available under the gpl v any restrictions to use by non-academics: none abbreviations chip-chip: chromatin immunoprecipitation followed by microarray; chip-seq: chromatin immunoprecipitation followed by sequencing; ctcf: ccctc binding factor; gabp: ga-binding protein; go: gene ontology; kde: kernel density estimation; macs: model-based analysis of chip-seq; nrsf: neuron-restrictive silencer transcription factor, also known as rest; pol ii: polymerase ii; rad21: homolog; taf1: tata box binding protein-associated factor 1; tf: transcription factor. competing interests the authors declare that they have no competing interests. authors’ contributions ll conceived the study, yl, and ll performed the analyses, and yl, ll, and dmu were all involved in study design, planning analyses, and the interpretation of results. all authors contributed to writing and revising the manuscript, and all approved the final manuscript for publication. supplementary material additional file 1: supplementary tables s1-19 cell lines contributing chip-seq data for each of the transcription factors. click here for file additional file 2 outlines of algorithms. t-kde. click here for file additional file 3: supplementary figure s1 proportion of tkde-declared versus kde-declared constitutive ctcf binding sites whose distance from nearest motif-based constitutive ctcf binding site on chromosomesare less than distance d plotted as a function of d for various bandwidths.separate curves for t-kde with bandwidth of bps and for the same density estimation algorithm without the binary range tree pre-processing. click here for file additional file 4 additional results comparing various methods. t-kde: a method for genome-wide identification of constitutive protein binding sites from multiple chip-seq data sets. click here for file acknowledgements we thank liang niu and weichun huang for discussion and xuting wang and grace kissling for critical reading of the manuscript. we thank the computational biology facility at niehs for computing time and support. this research was supported by intramural research program of the nih, national institute of environmental health sciences.
a protein may bind to its target dna sites constitutively, i.e., regardless of cell type. furthermore, t-kde can identify constitutive sites that are missed by a motif-based approach either because a bound site failed to reach the motif significance cutoff or because the peak sequence scanned was too short. to locate constitutive binding sites for a protein using chip-seq data for that protein from multiple cell lines, we developed a method, t-kde, which combines a binary range tree with a kernel density estimator. binding patternchip-seqkernel density estimationbinary range treemode-findingconstitutive sitectcf code t-kde is an efficient and effective method to predict constitutive protein binding sites using chip-seq peaks from multiple cell lines. using ctcf chip-seq datasets, we showed that the number of constitutive sites identified by t-kde is robust to the choice of tuning parameter and that t-kde identifies binding site locations more accurately than a binning approach. genome-wide discovery of constitutive binding sites requires robust and efficient computational methods to integrate results from numerous binding experiments.
 a protein may bind to its target dna sites constitutively, i.e., regardless of cell type. intuitively, constitutive binding sites should be biologically functional. a prerequisite for understanding their functional relevance is knowing all their locations for a protein of interest. genome-wide discovery of constitutive binding sites requires robust and efficient computational methods to integrate results from numerous binding experiments. such methods are lacking, however.  to locate constitutive binding sites for a protein using chip-seq data for that protein from multiple cell lines, we developed a method, t-kde, which combines a binary range tree with a kernel density estimator. using ctcf chip-seq datasets, we showed that the number of constitutive sites identified by t-kde is robust to the choice of tuning parameter and that t-kde identifies binding site locations more accurately than a binning approach. furthermore, t-kde can identify constitutive sites that are missed by a motif-based approach either because a bound site failed to reach the motif significance cutoff or because the peak sequence scanned was too short. by studying sites declared constitutive by t-kde but not by the motif-based approach, we discovered two new ctcf motif variants. using encode data on transcription factors in cell lines, we identified constitutive binding sites for each tf and provide evidence that, for some tfs, they may be biologically meaningful.  t-kde is an efficient and effective method to predict constitutive protein binding sites using chip-seq peaks from multiple cell lines. besides constitutive binding sites for a given protein, t-kde can identify genomic “hot spots” where several different proteins bind and, conversely, cell-type-specific sites bound by a given protein. binding patternchip-seqkernel density estimationbinary range treemode-findingconstitutive sitectcf code 