BACKGROUND
to reconstruct and model transcriptional regulatory networks  we need to know the genome-wide binding sites of transcription factors
<cit> .
chromatin immunoprecipitation followed by microarray
<cit>  or sequencing
<cit>  has been extensively used to study the in vivo binding locations of individual transcription factors and cofactors in a wide range of species and tissues
<cit> .
despite their popularity, such methods can only study a single tf in a single cell type, under a specific condition, in each experiment.
thus, it is difficult to use these methods to obtain a comprehensive understanding of the complicated mammalian trns.
these networks can involve hundreds or thousands of tfs whose activities change across different tissues and conditions.
using computational methods to integrate other genomic resources in order to predict tissue-specific transcription factor binding is therefore an important research challenge.
several methods have been developed to use in vitro data characterizing tf binding specificities to identify tf binding sites across the genome.
specifically, data from universal protein binding microarray
<cit>  is often used for such analysis.
pbm is capable of analyzing the interaction of a sequence-specific tf with tens of thousands of short dna sequences  in a single experiment, and thus provides a highly detailed picture of tf-dna interactions.
it has been successfully applied to reveal the binding profiles of hundreds of tfs in yeast
<cit> , worm
<cit> , mouse
<cit>  and arabidopsis
<cit> .
some of the proposed methods for using pbm data represent tf binding preference by position weight matrices
<cit> .
however, pwms, although popular due to their simplicity, assume independence between positions, an assumption which may not hold in many cases
<cit> .
in contrast, more sophisticated models  may better represent the full binding profiles of tfs, without loss of information from using pwm.
for instance, it has been suggested that many tfs have more than one binding preference
<cit>  and these are easier to represent using k-mers.
while in vitro data provides important information regarding binding specificities, such data is context independent.
actual binding is highly dependent on tissue-specific conditions including chromatin accessibility, the presence of co-factors, etc
<cit> .
recently, a number of studies have reported that epigenetic information including certain histone modifications and hypersensitivity to dnase i cleavage correlate with tf binding in vivo <cit> .
moreover, functional tfbss tend to be under stronger negative selection, leaving a "phylogenetic footprint" in the genomic sequences.
several methods for predicting in vivo tf binding sites have attempted to combine such information with pwms to predict global binding preferences
<cit> .
however, as mentioned above, pwm may not be the best representation of tf binding.
as we show, by using a model that retains the dependence between positions in the motif we can improve upon methods that integrate epigenetic and pwm data.
in addition, none of these methods have so far been applied to elucidate the complete set of targets for tfs across a large number of tissues.
to predict accurate tissue-specific tfbs, we integrate multiple types of genomic data.
the first part of our model is a biophysically-motivated k-mer based method for analyzing pbm data, which allows secondary binding profiles and nucleotide dependencies in different positions of the tf binding sites.
next, we develop a new method, pipes , to combine the results from the pbm model with dnase i hypersensitivity  data and evolutionary conservation data to predict tissue-specific tfbs in vivo.
we demonstrate that such an integrative model significantly boosts context specific prediction results compared with using pbm data alone.
we also show that pipes improves upon other methods developed for integrating data to predict tfbs
<cit> , in some cases significantly so.
finally, we created a resource for tissue-specific trns using pbm data for  <dig> mouse tfs from uniprobe
<cit>  and dnase i hypersensitivity data for  <dig> mouse tissue/cell types from the mouse encode project
<cit> .
we predict the activities of tfs across different tissues, and, as we show, many of these predictions agree with current knowledge regarding tissue-specific roles of tfs.
our tissue specific activity predictions are also supported by global analysis of tf expression data.
the comprehensive resource of tf binding sites we built thus provides a reference map for understanding complex gene expression patterns.
RESULTS
an overview of our pipes method is shown in figure
<dig>  our model has two components: the left part of the figure shows our model for the pbm data, and the right part our model of epigenetics and conservation data.
starting with raw fluorescent intensities measured by pbm, we first infer binding probabilities to each individual k-mer with a biophysically-motivated model .
this information, based on pbm alone, can be used to score a sequence for potential tfbs .
next, we use tissue specific dnase i hypersensitivity data to determine chromatin accessibility , and combine such information with sequence conservation and the pbm derived scores to predict in vivo binding sites .
k-mer based pbm analysis can accurately infer tf binding specificities
while a number of methods have been suggested to use pbm data for predicting tf binding sites , we decided to extend k-mer based methods using a biophysically-motivated model .
k-mer based methods were shown to achieve the best performance among several techniques for the analysis of pbm data
<cit> .
such methods allow an intuitive representation of potential alternative binding motifs and can account for dependency among positions in a motif.
we use lasso regression with positive constraints to learn model parameters that represent binding probabilities to individual k-mers, where k is determined as part of the learning procedure.
this results in a sparse model with relatively few k-mers having nonzero binding probabilities.
on average, the number of k-mers with non-zero probabilities is  <dig>  ±  <dig>  across  <dig> mouse tfs with pbm data available.
the model combines the benefits of recent pwm-based biophysical methods  with the ability of pbms to capture dependencies between positions in a given motif .
we illustrate the results of our pbm model using four tfs including sox <dig>  esrra, klf <dig> and pou2f <dig>  figure
<dig> presents the pwms derived from the pbm data for these tfs by the seed-and-wobble algorithm
<cit>  , pwms in transfac
<cit>  for the corresponding tfs when available, and all the k-mers estimated by our model that have binding probabilities above  <dig> .
for s&w pwms, when a secondary binding preference was derived
<cit> , both the primary and secondary pwms are shown.
as can be seen, our k-mer model does well for this data.
for sox <dig>  the learned k-mers match well with both the primary and secondary s&w pwms .
for esrra and klf <dig>  k-mers matching the consensus sequences of the primary s&w pwms, respectively, are predicted to have high binding probabilities , whereas k-mers matching the secondary s&w pwms have lower predicted binding probabilities ranging from  <dig>  to  <dig>  .
in the case of pou2f <dig>  none of the inferred top k-mers match the s&w pwm .
however, many of these top k-mers closely match the consensus sequence of the transfac motif for pou2f <dig> derived from literature evidence .
overall, these results support the use of a biophysically-motivated model: the binding probabilities of k-mers are largely consistent with the results from independent methods and known motifs from transfac.
to test our pbm model, and as a baseline, we next used the inferred binding probabilities to predict in vivo tf binding.
we collected  <dig> published mouse chip-seq datasets for which the pbm data for the same tf or for a tf with a similar dna-binding domain is available .
from each chip-seq dataset, the top  <dig> peaks with highest enrichment are extracted, and the 600bp genomic regions centered on the reported peaks are used as the positive sequences bound by the tf.
then, 600bp sequences that  are upstream of and  300bp apart from each positive sequence, and  do not overlap with any other positive sequences, are used as negative sequences.
we also explored two alternative options for constructing negative sequences, including using size-matched random promoter sequences and randomly generated sequences.
the choice of negative sequence sets makes very little difference on the auc values, so we will only report the results based on the first negative set here .
we compared our pbm model with seven other methods that predict affinities of tf binding to given sequences, and the different methods are evaluated using areas under the roc curve  as a measure of their abilities to correctly classify the two sets of sequences .
the results indicate that the performance of our pbm method is at least comparable, and in some cases better than, previous methods.
for  <dig> of the  <dig> tfs we tested , our method improved over all other methods.
in all other cases, the auc of our method ranks within the top  <dig> .
notably, in such cases, none of the other methods consistently achieves the best auc.
the pwm-based and e-score based methods tend to work well for some cases , but for others their performance is not as good.
overall, our pbm model has the highest average auc over the  <dig> tfs tested .
to further assess the benefits of using pbm to derive tf binding specificities, we use two other collections of pwms for auc evaluation.
the homer pwms
<cit>  were derived from chip-seq datasets, while the jaspar pwms were from multiple sources .
overall, the results from the jaspar pwms are very similar to those obtained from the seed and wobble pwms, and both are weaker than our method.
the use of homer pwms lead to better overall performance.
however, given that homer trains the pwms from chip-seq data and the same datasets may be used for evaluation, this is clearly not a fair comparison.
we did notice that for some tfs, our method outperforms homer pwms .
the full details are shown in supplementary results in additional file
<dig> and additional file
<dig>
integrated model of pbm and dnase i hypersensitivity data significantly improves tfbs prediction accuracy
pbm data, although powerful, only measures in vitro binding.
therefore, even when using sophisticated methods, the ability to predict in vivo binding based on pbm data is limited.
dnase i hypersensitive  sites are regions of chromatin that are very sensitive to dnase i cleavage
<cit> , and previous studies have shown that such hypersensitivity correlates with tf binding
<cit> .
to better predict tissue-specific in vivo binding sites, we developed pipes, a probabilistic graphical model for integrating dnase i hs data with pbm data.
for windows containing a 36bp genomic region , we assume the chromatin of the site could exist in two states: open or closed, and that only in the open state the chromatin is accessible to binding by a tf.
we infer the chromatin state by using a mixture model for the dnase hs data: the open state should be associated with higher tag densities from the dnase data, and the closed state with lower densities.
the in vivo occupancy of a site is then estimated as the probability of binding in vitro estimated using result from the pbm model multiplied by the probability that the site is in an open state inferred from the dnase hs data .
figure
<dig> presents the aucs from applying pipes to predict in vivo tf binding in the corresponding tissues for the same  <dig> tfs studied in the previous section.
compared with using pbm data alone , the incorporation of dnase i hs data in the corresponding tissues  improves performance for  <dig> of the  <dig> tfs.
overall, the improvement of auc from adding dnase hs data across all tfs is statistically significant .
the biggest improvements are seen for tfs for which the results when using only pbm data are relatively poor.
for example, when predicting srf binding sites in heart, even the best methods analyzed above achieve an auc only slightly better than random .
by integrating pbm and dnase data, the performance of our method is improved by 46% from  <dig>  to  <dig>  .
similar improvement is also observed for oct <dig> .
as another baseline, we evaluated the aucs from using dnase hs data alone.
somewhat surprisingly, this feature alone seems quite discriminative : the mean auc is  <dig> .
nevertheless, in  <dig> out of  <dig> cases, our model using both dnase and pbm data improves these baseline results and its mean auc is also higher at  <dig> .
we also point out that in practice it is not appropriate to use the hs data alone to predict binding sites for a tf as the predictions would not be specific to the tf of interest.
<dig> for more details.
in addition to dnase i hypersensitivity data, bona fide tf binding sites are usually under evolutionary pressure and therefore more conserved
<cit> .
we thus further extended pipes to incorporate phastcons scores
<cit>  for each site .
performance of the full model that incorporates phastcons information is shown in figure
<dig> .
as can be seen, while in some cases adding the conservation information very slightly improves performance , overall using conservation data does not lead to a significant improvement in prediction accuracy.
when the dnase hs data is not available, using phastcons in addition to the pbm data provides a slight improvement of the aucs from using pbm alone .
the average auc is increased by 1%, and for some tfs, the improvement can be quite significant .
finally, we compare pipes with recent methods proposed for integrating dnase and motif information to predict tfbs.
the first method we compare against, termed 'intersection", was used by neph et al.
<cit> .
this method intersects sites that have high-scoring pwm matches for a tf with dnase hs sites to predict in vivo binding .
we also compare pipes with centipede
<cit> .
centipede uses a probabilistic model to integrate the prior information of putative sites, such as sequence conservation and matches to pwms, with the epigenetic data to predict binding sites.
while probabilistic, centipede does, however, rely on a stringent cutoff for pwm match scores to achieve low false positive rates.
the results are presented in figure
<dig>  as can be seen, the intersection method leads to aucs that are significantly lower than the ones obtained by our method for all  <dig> tfs , indicating that strict cutoffs  may lead to a high rate of false negatives.
similarly, using the default settings for centipede led to auc scores that are much lower than ours .
to further explore this, we varied the setting of centipede, including using different pwms and a range of cutoffs for defining putative binding sites, but the results remained the same .
the difference in aucs is highly significant: p =  <dig>  · 10- <dig>  using the best combination of pwm and cutoff for centipede.
these results indicate that our pipes model, which relies on k-mer based representation and avoids strict cutoffs, can improve in vivo predictions of tfbs.
we performed additional analysis to the integration model.
for these, we replace the binding probabilities predicted from the k-mer model, with those predicted when using pwms.
three different versions of pwms were used: the seed-and-wobble pwms and rap pwms learned from the same pbm data, and the jaspar pwms.
in all cases, the aucs are substantially higher than the ones from the intersection method and centipede .
these results indicate that the probabilistic integration step alone is enough to improve upon prior methods .
the usefulness of the k-mer based analysis provides additional advantage, as independently demonstrated in the earlier section.
combining pbm and dnase data enables the prediction of tissue-specific tf activities
the recently released mouse encode project data provides dnase i hypersensitivity data for more than  <dig> mouse tissue/cell types .
we set out to combine the pbm data for  <dig> mouse tfs in uniprobe with such dnase data to predict tissue-specific tf targets and determine tissue-specific tf activities .
identifying tfs that are highly active in specific tissues is useful for determining the function of such tfs, and serves as an initial step for reconstructing the tissue-specific transcriptional regulatory networks.
we predict how likely a tf is functional in any given tissue/cell type with an activity score for each tf-tissue pair.
our hypothesis is that if the tf is active in a tissue, it will bind a number of target sequences, thus the putative tf binding sites will be overrepresented in the dnase hs regions .
a higher activity score indicates that the tf is more active in the corresponding tissue .
we use a binomial test to assess the statistical significance of the activity scores .
the complete results, including activity scores and p-values, are provided in additional file
<dig>  in figure
4a we illustrate these results by focusing on the activity scores calculated for  <dig> tfs  across  <dig> representative tissue/cell types.
gata <dig> is known to function in mouse fetal liver haematopoiesis
<cit> , and its expression had also been observed in leukemia cells
<cit> .
our results are in good agreement with the prior knowledge regarding gata3’s activity: the top two tissues predicted for gata <dig> are e <dig>  liver cells and the adult leukemia cell line.
similarly the top tissue for pou6f <dig> is e <dig>  whole brain, in agreement with its known role in brain development
<cit> .
crx is an important tf for regulating photoreceptor genes in retina
<cit> , and our method correctly determined that its activity score in that tissue is the highest.
finally, hnf4a is a well known master regulator of liver- and kidney-specific genes
<cit> , as correctly predicted by our method.
while we only show  <dig> tissues, for all four tfs the correct tissues shown in figure
4a have the highest scores among all  <dig> tissues we tested .
to more globally validate these tissue-specific tf activities, we compared the correlation between our predicted tf activity scores and mrna levels for the same tfs in the corresponding tissue .
eight tissues and  <dig> tfs that are common to both datasets are used .
even though the two types of data  measure completely different aspects of cellular activity, we observe a pearson correlation coefficient of  <dig> , which is highly statistically significant .
since many tfs are only post-transcriptionally regulated, such a significant correlation provides strong support to the predictions computed by our method.
existing literature strongly supports predicted tf activities in several tissues
to further validate our predictions and investigate their potentials to lead to new biological insights, we took a closer look at the tfs predicted to be active in the adult liver tissue.
the top five such predictions are shown in table
1a .
besides hnf4a discussed above, rara, nr2f <dig>  rxra and tcf <dig> are all known to either regulate liver-specific genes or are involved in maintaining liver metabolism and homeostasis .
the 7th ranked factor tcf7l <dig>  was linked to type  <dig> diabetes risk in previous studies using snp data
<cit> , but the mechanism for its involvement was unclear.
our result indicates that it may have a regulatory role in liver metabolism.
indeed, a very recent study confirms its role in regulating key liver-specific metabolic genes
<cit> .
our result also assign a high liver activity score to cutl <dig> .
cutl <dig> was a known transcriptional repressor of terminal differentiation genes in several cell lineages including hepatocyte
<cit> .
recently, cutl <dig> was identified as target of the liver-specific microrna mir <dig> and a central mediator of the effects caused by the deregulation of mir <dig> in hepatocellular carcinoma
<cit> .
further down the list, foxa <dig>  is known to regulate lipid metabolism and ketogenesis related genes in liver
<cit> , and lef <dig>  is a prognostic biomarker for liver metastasis in colorectal cancer.
other tfs ranked within the top  <dig> for liver include tcf <dig> and tcf <dig>  members of the t-cell factor  family that are critical for hepatocyte metabolism and function
<cit> ; bhlhb <dig>  which is involved in the regulation of lipogenesis in liver
<cit> ; and hmbox <dig>  whose expression levels was shown to be reduced in liver cancer compared with surrounding normal tissues
<cit> .
overall, our predicted set of liver regulators is comprehensive, spanning several different classes of liver related activities including glucose and lipid metabolism and cancer, and including both repressors and activators.
in addition, table
1b and c presents the top  <dig> predicted tfs for two more tissues .
as can be seen, for almost all of these tfs there is strong support for their tissue-specific activity in the predicted tissue.
discussion
a number of recent projects including encode
<cit> , modencode
<cit>  and the roadmap epigenomics project
<cit>  have generated large amounts of genomic data.
an important research goal is to translate these resources into accurate, tissue and condition-sensitive, molecular-level networks.
constructing tissue-specific maps of tf binding sites is a central part of this overall research effort.
using several different datasets and a novel computational strategy, pipes, we demonstrated that such high-quality computational predictions can be obtained.
we used pipes to compile a resource that includes comprehensive predictions for more than  <dig> tfs across  <dig> tissues.
a recent benchmark study that compared many methods for analyzing pbm data concluded that pwm-based methods work as well as other models for predicting tfbs
<cit> .
our results differ from these previous studies.
this could either be the result of the lasso based method we have used or the specific dataset we used for the comparison.
additional work is required to reach a definitive conclusion regarding the importance of independence assumption used by pwms when modeling tfbs.
here we focused on integrating a number of datasets for predicting tfbss.
for our biophysical approach, using a k-mer based method allowed us to capture both the dependency among positions within the binding site as well as multiple different motifs for a single tf.
such method worked well in classifying bound and unbound sequences from in vivo chip-seq data for many tfs.
moreover, integrating pbm data with chromatin accessibility from dnase i hs data greatly improves the accuracy of tf binding predictions.
several recent papers explored related ideas.
chromia
<cit>  used a hidden markov model to combine sequence-specific tf binding with histone modification data, but their predictions were based on pwm scoring and only focused on a dozen of tfs in mouse embryonic stem cells.
ernst et al.
<cit>  combined experimental data from a number of tissues to generate a single  tf-target prediction map.
however, that method has also relied on pwms and no tissue specific predictions were made.
centipede
<cit>  used unsupervised methods to integrate tf-dna interaction, epigenetic and evolutionary data, and is most similar to our efforts.
however, centipede relies on a footprint in the dnase data that tfs leave.
such dnase footprints are the actual locations where the tf binds and are therefore protected from dnase cleavage within the dnase hs site.
unfortunately, dnase footprint data is expensive to obtain  since it requires very high coverage when sequencing.
in addition, centipede uses a stringent cutoff  to define putative binding sites, and thus may lose significant information in relatively weak binding sites, which have been shown to be collectively important for tf binding
<cit> .
neph et al.
<cit>  also combined dnase footprints with pwms to predict tf-tf interactions  across a large number of human tissues, using a simple method to intersect motif matches with dhs sites.
such hard cutoffs may miss sites that score high  for both types of data which are found by our method.
in general, we find that binding sites for the tfs we looked at are only modestly conserved when compared with controls .
this is largely consistent with the recent findings that functional non-coding sequences evolve rather rapidly
<cit> .
as a result, adding phastcons in the integrated model does not lead to improvements in auc values.
nevertheless, there are a number of advantages for models that can incorporate sequence conservation.
first, when dnase data is not available, adding conservation in the model leads to slight improvement over models that only use pbm data .
second, the conservation of binding events can vary greatly among tissues.
for example, enhancers in brain are far more constrained than those in the heart
<cit> .
thus it is quite possible that sequence conservation would be more informative for other chip-seq studies.
the application of pipes to predict tissue specific tfbs led to results that agree well with existing knowledge regarding tf roles in specific tissues.
the overall results are significantly correlated with independent gene expression data measured for these tfs across tissues even though such expression data was not used at all in our analysis.
several extensions of our current work are possible.
recently jiang et al.
<cit>  reported the interesting phenomenon of sticky k-mers: these are k-mers that appear to bind to tfs with relatively high affinities in a large number of pbm experiments.
the sticky k-mers likely represent background noises in the pbm experiments and an interesting research direction is to expand our regression method to remove such noises.
in another recent study, ballare et al.
<cit>  reported that functional tfbss are not always associated with high chromatin accessibility, an assumption implicitly made by us and other related methods.
rather, nucleosomes may occupy tfbss at basal conditions, and are only remodeled or displaced upon change of cellular conditions .
despite this unexpected relationship between tfbs and chromatin states, the paper does report that such sites, while associated with high nucleosome occupancy before stimulation, often overlap with dnase hs sites.
it remains to be seen how common such cases are and what is the impact on methods such as ours that rely on dnase data to predict condition specific tf binding.
moreover, our integrative framework for utilizing additional information sources when predicting binding events on a genome wide scale could also be used for large scale comparison of different pwm methods and methods that use more complicated models to represent tf binding preferences
<cit> .
this requires a detailed study and is left for future work.
CONCLUSIONS
combining pbm and dnase data, we presented the first major effort to provide a systematic map of computationally predicted tissue-specific targets for hundreds of tfs across a large number of tissues in mouse.
we complied a resource that provides tf-target predictions for all  <dig> tfs studied across the  <dig> tissue/cell types .
we believe that such comprehensive resource would be useful for both biology and computation-oriented researchers studying gene regulation
<cit> .

 to reconstruct and model transcriptional regulatory networks we need to know the genome-wide binding sites of transcription factors. chromatin immunoprecipitation followed by microarray or sequencing has been extensively used to study the in vivo binding locations of individual transcription factors and cofactors in a wide range of species and tissues. despite their popularity, such methods can only study a single tf in a single cell type, under a specific condition, in each experiment. thus, it is difficult to use these methods to obtain a comprehensive understanding of the complicated mammalian trns. these networks can involve hundreds or thousands of tfs whose activities change across different tissues and conditions. using computational methods to integrate other genomic resources in order to predict tissue-specific transcription factor binding is therefore an important research challenge. several methods have been developed to use in vitro data characterizing tf binding specificities to identify tf binding sites across the genome. specifically, data from universal protein binding microarray is often used for such analysis. pbm is capable of analyzing the interaction of a sequence-specific tf with tens of thousands of short dna sequences in a single experiment, and thus provides a highly detailed picture of tf-dna interactions. it has been successfully applied to reveal the binding profiles of hundreds of tfs in yeast, worm, mouse and arabidopsis. some of the proposed methods for using pbm data represent tf binding preference by position weight matrices. however, pwms, although popular due to their simplicity, assume independence between positions, an assumption which may not hold in many cases. in contrast, more sophisticated models may better represent the full binding profiles of tfs, without loss of information from using pwm. for instance, it has been suggested that many tfs have more than one binding preference and these are easier to represent using k-mers. while in vitro data provides important information regarding binding specificities, such data is context independent. actual binding is highly dependent on tissue-specific conditions including chromatin accessibility, the presence of co-factors, etc. recently, a number of studies have reported that epigenetic information including certain histone modifications and hypersensitivity to dnase i cleavage correlate with tf binding in vivo. moreover, functional tfbss tend to be under stronger negative selection, leaving a "phylogenetic footprint" in the genomic sequences. several methods for predicting in vivo tf binding sites have attempted to combine such information with pwms to predict global binding preferences. however, as mentioned above, pwm may not be the best representation of tf binding. as we show, by using a model that retains the dependence between positions in the motif we can improve upon methods that integrate epigenetic and pwm data. in addition, none of these methods have so far been applied to elucidate the complete set of targets for tfs across a large number of tissues. to predict accurate tissue-specific tfbs, we integrate multiple types of genomic data. the first part of our model is a biophysically-motivated k-mer based method for analyzing pbm data, which allows secondary binding profiles and nucleotide dependencies in different positions of the tf binding sites. next, we develop a new method, pipes, to combine the results from the pbm model with dnase i hypersensitivity data and evolutionary conservation data to predict tissue-specific tfbs in vivo. we demonstrate that such an integrative model significantly boosts context specific prediction results compared with using pbm data alone. we also show that pipes improves upon other methods developed for integrating data to predict tfbs, in some cases significantly so. finally, we created a resource for tissue-specific trns using pbm data for mouse tfs from uniprobe and dnase i hypersensitivity data for mouse tissue/cell types from the mouse encode project. we predict the activities of tfs across different tissues, and, as we show, many of these predictions agree with current knowledge regarding tissue-specific roles of tfs. our tissue specific activity predictions are also supported by global analysis of tf expression data. the comprehensive resource of tf binding sites we built thus provides a reference map for understanding complex gene expression patterns.  an overview of our pipes method is shown in figure our model has two components: the left part of the figure shows our model for the pbm data, and the right part our model of epigenetics and conservation data. starting with raw fluorescent intensities measured by pbm, we first infer binding probabilities to each individual k-mer with a biophysically-motivated model. this information, based on pbm alone, can be used to score a sequence for potential tfbs. next, we use tissue specific dnase i hypersensitivity data to determine chromatin accessibility, and combine such information with sequence conservation and the pbm derived scores to predict in vivo binding sites. k-mer based pbm analysis can accurately infer tf binding specificities while a number of methods have been suggested to use pbm data for predicting tf binding sites, we decided to extend k-mer based methods using a biophysically-motivated model. k-mer based methods were shown to achieve the best performance among several techniques for the analysis of pbm data. such methods allow an intuitive representation of potential alternative binding motifs and can account for dependency among positions in a motif. we use lasso regression with positive constraints to learn model parameters that represent binding probabilities to individual k-mers, where k is determined as part of the learning procedure. this results in a sparse model with relatively few k-mers having nonzero binding probabilities. on average, the number of k-mers with non-zero probabilities is ± across mouse tfs with pbm data available. the model combines the benefits of recent pwm-based biophysical methods with the ability of pbms to capture dependencies between positions in a given motif. we illustrate the results of our pbm model using four tfs including sox esrra, klf and pou2f figure presents the pwms derived from the pbm data for these tfs by the seed-and-wobble algorithm, pwms in transfac for the corresponding tfs when available, and all the k-mers estimated by our model that have binding probabilities above. for s&w pwms, when a secondary binding preference was derived, both the primary and secondary pwms are shown. as can be seen, our k-mer model does well for this data. for sox the learned k-mers match well with both the primary and secondary s&w pwms. for esrra and klf k-mers matching the consensus sequences of the primary s&w pwms, respectively, are predicted to have high binding probabilities, whereas k-mers matching the secondary s&w pwms have lower predicted binding probabilities ranging from to. in the case of pou2f none of the inferred top k-mers match the s&w pwm. however, many of these top k-mers closely match the consensus sequence of the transfac motif for pou2f derived from literature evidence. overall, these results support the use of a biophysically-motivated model: the binding probabilities of k-mers are largely consistent with the results from independent methods and known motifs from transfac. to test our pbm model, and as a baseline, we next used the inferred binding probabilities to predict in vivo tf binding. we collected published mouse chip-seq datasets for which the pbm data for the same tf or for a tf with a similar dna-binding domain is available. from each chip-seq dataset, the top peaks with highest enrichment are extracted, and the 600bp genomic regions centered on the reported peaks are used as the positive sequences bound by the tf. then, 600bp sequences that are upstream of and 300bp apart from each positive sequence, and do not overlap with any other positive sequences, are used as negative sequences. we also explored two alternative options for constructing negative sequences, including using size-matched random promoter sequences and randomly generated sequences. the choice of negative sequence sets makes very little difference on the auc values, so we will only report the results based on the first negative set here. we compared our pbm model with seven other methods that predict affinities of tf binding to given sequences, and the different methods are evaluated using areas under the roc curve as a measure of their abilities to correctly classify the two sets of sequences. the results indicate that the performance of our pbm method is at least comparable, and in some cases better than, previous methods. for of the tfs we tested, our method improved over all other methods. in all other cases, the auc of our method ranks within the top. notably, in such cases, none of the other methods consistently achieves the best auc. the pwm-based and e-score based methods tend to work well for some cases, but for others their performance is not as good. overall, our pbm model has the highest average auc over the tfs tested. to further assess the benefits of using pbm to derive tf binding specificities, we use two other collections of pwms for auc evaluation. the homer pwms were derived from chip-seq datasets, while the jaspar pwms were from multiple sources. overall, the results from the jaspar pwms are very similar to those obtained from the seed and wobble pwms, and both are weaker than our method. the use of homer pwms lead to better overall performance. however, given that homer trains the pwms from chip-seq data and the same datasets may be used for evaluation, this is clearly not a fair comparison. we did notice that for some tfs, our method outperforms homer pwms. the full details are shown in supplementary results in additional file and additional file integrated model of pbm and dnase i hypersensitivity data significantly improves tfbs prediction accuracy pbm data, although powerful, only measures in vitro binding. therefore, even when using sophisticated methods, the ability to predict in vivo binding based on pbm data is limited. dnase i hypersensitive sites are regions of chromatin that are very sensitive to dnase i cleavage, and previous studies have shown that such hypersensitivity correlates with tf binding. to better predict tissue-specific in vivo binding sites, we developed pipes, a probabilistic graphical model for integrating dnase i hs data with pbm data. for windows containing a 36bp genomic region, we assume the chromatin of the site could exist in two states: open or closed, and that only in the open state the chromatin is accessible to binding by a tf. we infer the chromatin state by using a mixture model for the dnase hs data: the open state should be associated with higher tag densities from the dnase data, and the closed state with lower densities. the in vivo occupancy of a site is then estimated as the probability of binding in vitro estimated using result from the pbm model multiplied by the probability that the site is in an open state inferred from the dnase hs data. figure presents the aucs from applying pipes to predict in vivo tf binding in the corresponding tissues for the same tfs studied in the previous section. compared with using pbm data alone, the incorporation of dnase i hs data in the corresponding tissues improves performance for of the tfs. overall, the improvement of auc from adding dnase hs data across all tfs is statistically significant. the biggest improvements are seen for tfs for which the results when using only pbm data are relatively poor. for example, when predicting srf binding sites in heart, even the best methods analyzed above achieve an auc only slightly better than random. by integrating pbm and dnase data, the performance of our method is improved by 46% from to. similar improvement is also observed for oct. as another baseline, we evaluated the aucs from using dnase hs data alone. somewhat surprisingly, this feature alone seems quite discriminative: the mean auc is. nevertheless, in out of cases, our model using both dnase and pbm data improves these baseline results and its mean auc is also higher at. we also point out that in practice it is not appropriate to use the hs data alone to predict binding sites for a tf as the predictions would not be specific to the tf of interest. for more details. in addition to dnase i hypersensitivity data, bona fide tf binding sites are usually under evolutionary pressure and therefore more conserved. we thus further extended pipes to incorporate phastcons scores for each site. performance of the full model that incorporates phastcons information is shown in figure. as can be seen, while in some cases adding the conservation information very slightly improves performance, overall using conservation data does not lead to a significant improvement in prediction accuracy. when the dnase hs data is not available, using phastcons in addition to the pbm data provides a slight improvement of the aucs from using pbm alone. the average auc is increased by 1%, and for some tfs, the improvement can be quite significant. finally, we compare pipes with recent methods proposed for integrating dnase and motif information to predict tfbs. the first method we compare against, termed'intersection", was used by neph et al.. this method intersects sites that have high-scoring pwm matches for a tf with dnase hs sites to predict in vivo binding. we also compare pipes with centipede. centipede uses a probabilistic model to integrate the prior information of putative sites, such as sequence conservation and matches to pwms, with the epigenetic data to predict binding sites. while probabilistic, centipede does, however, rely on a stringent cutoff for pwm match scores to achieve low false positive rates. the results are presented in figure as can be seen, the intersection method leads to aucs that are significantly lower than the ones obtained by our method for all tfs, indicating that strict cutoffs may lead to a high rate of false negatives. similarly, using the default settings for centipede led to auc scores that are much lower than ours. to further explore this, we varied the setting of centipede, including using different pwms and a range of cutoffs for defining putative binding sites, but the results remained the same. the difference in aucs is highly significant: p =   · 10- using the best combination of pwm and cutoff for centipede. these results indicate that our pipes model, which relies on k-mer based representation and avoids strict cutoffs, can improve in vivo predictions of tfbs. we performed additional analysis to the integration model. for these, we replace the binding probabilities predicted from the k-mer model, with those predicted when using pwms. three different versions of pwms were used: the seed-and-wobble pwms and rap pwms learned from the same pbm data, and the jaspar pwms. in all cases, the aucs are substantially higher than the ones from the intersection method and centipede. these results indicate that the probabilistic integration step alone is enough to improve upon prior methods. the usefulness of the k-mer based analysis provides additional advantage, as independently demonstrated in the earlier section. combining pbm and dnase data enables the prediction of tissue-specific tf activities the recently released mouse encode project data provides dnase i hypersensitivity data for more than mouse tissue/cell types. we set out to combine the pbm data for mouse tfs in uniprobe with such dnase data to predict tissue-specific tf targets and determine tissue-specific tf activities. identifying tfs that are highly active in specific tissues is useful for determining the function of such tfs, and serves as an initial step for reconstructing the tissue-specific transcriptional regulatory networks. we predict how likely a tf is functional in any given tissue/cell type with an activity score for each tf-tissue pair. our hypothesis is that if the tf is active in a tissue, it will bind a number of target sequences, thus the putative tf binding sites will be overrepresented in the dnase hs regions. a higher activity score indicates that the tf is more active in the corresponding tissue. we use a binomial test to assess the statistical significance of the activity scores. the complete results, including activity scores and p-values, are provided in additional file in figure 4a we illustrate these results by focusing on the activity scores calculated for tfs across representative tissue/cell types. gata is known to function in mouse fetal liver haematopoiesis, and its expression had also been observed in leukemia cells. our results are in good agreement with the prior knowledge regarding gata3’s activity: the top two tissues predicted for gata are e liver cells and the adult leukemia cell line. similarly the top tissue for pou6f is e whole brain, in agreement with its known role in brain development. crx is an important tf for regulating photoreceptor genes in retina, and our method correctly determined that its activity score in that tissue is the highest. finally, hnf4a is a well known master regulator of liver- and kidney-specific genes, as correctly predicted by our method. while we only show tissues, for all four tfs the correct tissues shown in figure 4a have the highest scores among all tissues we tested. to more globally validate these tissue-specific tf activities, we compared the correlation between our predicted tf activity scores and mrna levels for the same tfs in the corresponding tissue. eight tissues and tfs that are common to both datasets are used. even though the two types of data measure completely different aspects of cellular activity, we observe a pearson correlation coefficient of, which is highly statistically significant. since many tfs are only post-transcriptionally regulated, such a significant correlation provides strong support to the predictions computed by our method. existing literature strongly supports predicted tf activities in several tissues to further validate our predictions and investigate their potentials to lead to new biological insights, we took a closer look at the tfs predicted to be active in the adult liver tissue. the top five such predictions are shown in table 1a. besides hnf4a discussed above, rara, nr2f rxra and tcf are all known to either regulate liver-specific genes or are involved in maintaining liver metabolism and homeostasis. the 7th ranked factor tcf7l was linked to type diabetes risk in previous studies using snp data, but the mechanism for its involvement was unclear. our result indicates that it may have a regulatory role in liver metabolism. indeed, a very recent study confirms its role in regulating key liver-specific metabolic genes. our result also assign a high liver activity score to cutl. cutl was a known transcriptional repressor of terminal differentiation genes in several cell lineages including hepatocyte. recently, cutl was identified as target of the liver-specific microrna mir and a central mediator of the effects caused by the deregulation of mir in hepatocellular carcinoma. further down the list, foxa is known to regulate lipid metabolism and ketogenesis related genes in liver, and lef is a prognostic biomarker for liver metastasis in colorectal cancer. other tfs ranked within the top for liver include tcf and tcf members of the t-cell factor family that are critical for hepatocyte metabolism and function; bhlhb which is involved in the regulation of lipogenesis in liver; and hmbox whose expression levels was shown to be reduced in liver cancer compared with surrounding normal tissues. overall, our predicted set of liver regulators is comprehensive, spanning several different classes of liver related activities including glucose and lipid metabolism and cancer, and including both repressors and activators. in addition, table 1b and c presents the top predicted tfs for two more tissues. as can be seen, for almost all of these tfs there is strong support for their tissue-specific activity in the predicted tissue. discussion a number of recent projects including encode, modencode and the roadmap epigenomics project have generated large amounts of genomic data. an important research goal is to translate these resources into accurate, tissue and condition-sensitive, molecular-level networks. constructing tissue-specific maps of tf binding sites is a central part of this overall research effort. using several different datasets and a novel computational strategy, pipes, we demonstrated that such high-quality computational predictions can be obtained. we used pipes to compile a resource that includes comprehensive predictions for more than tfs across tissues. a recent benchmark study that compared many methods for analyzing pbm data concluded that pwm-based methods work as well as other models for predicting tfbs. our results differ from these previous studies. this could either be the result of the lasso based method we have used or the specific dataset we used for the comparison. additional work is required to reach a definitive conclusion regarding the importance of independence assumption used by pwms when modeling tfbs. here we focused on integrating a number of datasets for predicting tfbss. for our biophysical approach, using a k-mer based method allowed us to capture both the dependency among positions within the binding site as well as multiple different motifs for a single tf. such method worked well in classifying bound and unbound sequences from in vivo chip-seq data for many tfs. moreover, integrating pbm data with chromatin accessibility from dnase i hs data greatly improves the accuracy of tf binding predictions. several recent papers explored related ideas. chromia used a hidden markov model to combine sequence-specific tf binding with histone modification data, but their predictions were based on pwm scoring and only focused on a dozen of tfs in mouse embryonic stem cells. ernst et al. combined experimental data from a number of tissues to generate a single tf-target prediction map. however, that method has also relied on pwms and no tissue specific predictions were made. centipede used unsupervised methods to integrate tf-dna interaction, epigenetic and evolutionary data, and is most similar to our efforts. however, centipede relies on a footprint in the dnase data that tfs leave. such dnase footprints are the actual locations where the tf binds and are therefore protected from dnase cleavage within the dnase hs site. unfortunately, dnase footprint data is expensive to obtain since it requires very high coverage when sequencing. in addition, centipede uses a stringent cutoff to define putative binding sites, and thus may lose significant information in relatively weak binding sites, which have been shown to be collectively important for tf binding. neph et al. also combined dnase footprints with pwms to predict tf-tf interactions across a large number of human tissues, using a simple method to intersect motif matches with dhs sites. such hard cutoffs may miss sites that score high for both types of data which are found by our method. in general, we find that binding sites for the tfs we looked at are only modestly conserved when compared with controls. this is largely consistent with the recent findings that functional non-coding sequences evolve rather rapidly. as a result, adding phastcons in the integrated model does not lead to improvements in auc values. nevertheless, there are a number of advantages for models that can incorporate sequence conservation. first, when dnase data is not available, adding conservation in the model leads to slight improvement over models that only use pbm data. second, the conservation of binding events can vary greatly among tissues. for example, enhancers in brain are far more constrained than those in the heart. thus it is quite possible that sequence conservation would be more informative for other chip-seq studies. the application of pipes to predict tissue specific tfbs led to results that agree well with existing knowledge regarding tf roles in specific tissues. the overall results are significantly correlated with independent gene expression data measured for these tfs across tissues even though such expression data was not used at all in our analysis. several extensions of our current work are possible. recently jiang et al. reported the interesting phenomenon of sticky k-mers: these are k-mers that appear to bind to tfs with relatively high affinities in a large number of pbm experiments. the sticky k-mers likely represent background noises in the pbm experiments and an interesting research direction is to expand our regression method to remove such noises. in another recent study, ballare et al. reported that functional tfbss are not always associated with high chromatin accessibility, an assumption implicitly made by us and other related methods. rather, nucleosomes may occupy tfbss at basal conditions, and are only remodeled or displaced upon change of cellular conditions. despite this unexpected relationship between tfbs and chromatin states, the paper does report that such sites, while associated with high nucleosome occupancy before stimulation, often overlap with dnase hs sites. it remains to be seen how common such cases are and what is the impact on methods such as ours that rely on dnase data to predict condition specific tf binding. moreover, our integrative framework for utilizing additional information sources when predicting binding events on a genome wide scale could also be used for large scale comparison of different pwm methods and methods that use more complicated models to represent tf binding preferences. this requires a detailed study and is left for future work.  combining pbm and dnase data, we presented the first major effort to provide a systematic map of computationally predicted tissue-specific targets for hundreds of tfs across a large number of tissues in mouse. we complied a resource that provides tf-target predictions for all tfs studied across the tissue/cell types. we believe that such comprehensive resource would be useful for both biology and computation-oriented researchers studying gene regulation. 
we demonstrate that pipes improves over existing methods on distinguishing between in vivo bound and unbound sequences using chip-seq data for mouse tfs. we provide a systematic map of computationally predicted tissue-specific binding targets for mouse tfs across tissue/cell types. such comprehensive resource is useful for researchers studying gene regulation. in addition, our predictions are in good agreement with current knowledge of tissue-specific tf regulation. pipes integrates in vitro protein binding microarrays, sequence conservation and tissue-specific epigenetic information. experimental methods to study in vivo binding, including chip-chip and chip-seq, can only study one transcription factor in a single cell type and under a specific condition in each experiment, and therefore cannot scale to determine the full set of regulatory interactions in mammalian transcriptional regulatory networks.
 studies of gene regulation often utilize genome-wide predictions of transcription factor binding sites. most existing prediction methods are based on sequence information alone, ignoring biological contexts such as developmental stages and tissue types. experimental methods to study in vivo binding, including chip-chip and chip-seq, can only study one transcription factor in a single cell type and under a specific condition in each experiment, and therefore cannot scale to determine the full set of regulatory interactions in mammalian transcriptional regulatory networks.  we developed a new computational approach, pipes, for predicting tissue-specific tf binding. pipes integrates in vitro protein binding microarrays, sequence conservation and tissue-specific epigenetic information. we demonstrate that pipes improves over existing methods on distinguishing between in vivo bound and unbound sequences using chip-seq data for mouse tfs. in addition, our predictions are in good agreement with current knowledge of tissue-specific tf regulation.  we provide a systematic map of computationally predicted tissue-specific binding targets for mouse tfs across tissue/cell types. such comprehensive resource is useful for researchers studying gene regulation. 