BACKGROUND
understanding regulatory circuits controlling gene expression is one of the fundamental problems in modern biology.
gene expression is controlled at many different levels but control of transcription is one of the main steps of regulation.
one of the best understood control mechanisms is the binding of transcription factors  to the regulatory sites on dna in a sequence-specific manner, which affects transcription initiation  <cit> .
the important problem of locating the binding sites for specific tfs, and thus identifying the genes they regulate, has attracted much attention from the bioinformatics community  <cit> .
different methods have been employed for abstracting patterns or "motifs" from the sequences that bind particular tfs leading to predictions of likely binding sites in the genome of the organism under study.
factors regulating multiple genes often have binding motifs low in information content  <cit> , making the task of prediction harder.
examples of such highly pleiotropic proteins range from global regulators in prokaryotes  to hox proteins  <cit> , important in metazoan development.
experimental approaches to locating binding sites on dna  <cit> , have uncovered numerous binding sites for various factors.
however, looking at the databases devoted to such regulatory sites, like dpinteract  <cit>  and regulondb  <cit>  for e. coli, scpd for yeast  <cit>  and transfac for many higher eukaryotic organisms  <cit> , it is obvious that, for most pleiotropic tfs targeting a large number  of genes, the number of known sites is still a small fraction of all the functional sites.
a high-throughput version of the chromatin immunoprecipitation method, commonly known as the "chip on chip", has been introduced recently  <cit> .
in principle, this method locates binding sites genome-wide.
however, the resolution is limited to several hundred bases and requires further bioinformatic analysis  <cit> .
an alternative approach would be to find the dna binding specificity of a tf by an in vitro method and then use the binding motif to search the genome for putative sites.
one of these methods is selex  <cit> , which is often used to find the strongest binding sites  from a library consisting of randomly generated oligonucleotides.
however, a tf can often function at binding sites that are far weaker than the consensus.
therefore, to characterize the binding preferences of a tf, we need to identify many of these potential weak binding sites and to estimate the parameters describing the statistical distribution of those sequences.
the appropriate modification of the selex procedure needed to achieve this goal is based on the selex-sage procedure  <cit> .
analysis of the conditions under which we get a significant number of intermediate strength sites was performed in  <cit> .
we will use this procedure on the pleiotropic e. coli factor cap.
an alternative to this technology would have been to use dna chips for protein binding  <cit> .
currently, for transcription factors with long binding sites , it is common practice to use genomic sequences rather than random libraries in dna chips.
this has its advantages but also might lead to uncertainties regarding the genomic background model in the final statistical analysis.
to abstract a motif from the sequences found by the modified selex process, we need a computational method: a supervised algorithm, trained on a set of binding sites identified directly by experimental measurements  <cit> .
we will compare different supervised methods for extraction of parameters and use cap targets as a benchmark.
the widely used bioinformatic tool for quantitatively describing such motifs is the weight matrix method  <cit> .
setting the threshold correctly is essential for the quality of predictions .
however, optimization of the threshold is a non-trivial problem, resolving which is one of the goals of this study.
we have shown  <cit>  that using the physically correct expression for binding probability, with saturation effects built in, leads to a more accurate estimate for the binding energy and provides a practically useful solution to the problem of classifier threshold choice.
the resulting method, quadratic programming method of energy matrix estimation or qpmeme  <cit> , turns out to be a one-class support vector machine  <cit> .
in this paper we do the following:
perform high throughput, low stringency selex experiments.
analyze selex experiments and extract parameters for models of sequence-dependent tf/dna interactions by using qpmeme and its extensions.
verify in vitro affinity predicted by models built on selex data by electrophoretic mobility shift assay.
identify cap targets in e. coli.
compare various predicted sites from different methods to see their functional conservation among gram-negative facultative anaerobes.
RESULTS
sage selex study of cap
we followed the selex-sage procedure  <cit> .
after extracting a large number of sequences, we passed these sequences through a quality filter to find sequences that were likely to have been bound by a cap dimer in the random region of the selex pool.
we were left with  <dig> sequences at the end.
analysis of sequences from sage selex study
the maximum-likelihood method for distributions with sharp cutoffs was described in  <cit> .
in particular circumstances, the parameter estimation method becomes a support vector machine, and the resulting algorithm, qpmeme, can be used to determine the binding energy of the protein to the sequence s, as e = ∑iα εiα siα .
the variable siα is defined as follows: if the ith base is α in sequence s, then siα =  <dig> and siα =  <dig> otherwise.
the parameter εiα is the contribution to the binding free energy from base α being at position i. these εiα parameters are chosen to minimize the variance of e over the background distribution of sequences s, subject to the constraints e) ≤ - <dig> for the set of example binding sequences s, j =  <dig>  ..., n. sequences satisfying e) ≤ - <dig> are then declared to be binding sites.
in practice, the base frequencies are taken to be independent and the probability of finding the base α is taken to be pα.
the quantity to be minimized is given by ∑iαpαεiα <dig>  subject to the constraints ∑αpαεiα= <dig>  for each i.
the qpmeme algorithm was used to produce an energy matrix, , using the set of  <dig> known cap binding sites from the dpinteract database.
an energy matrix was also constructed using the binding sites identified by the selex procedure described below.
weight matrices were also constructed for both the known sites and the selex sites using the formula wiα = log , where wiα is the i, α component of the weight matrix and fiα=∑jsiα is the frequency of the ith base being α. the background probabilities for g and c are taken to be the same.
the same applies to the background probabilities for a and t. for both algorithms, the background gc content was taken to be  <dig>   when constructing energy matrices from the known sites, and was taken to be  <dig>  when constructing energy matrices from selex sites.
for the weight matrices/energy matrices used in the study, see additional file  <dig>  additional file  <dig>  additional file  <dig> and additional file  <dig>
prediction of in vitro affinity of binding sites
the energy matrices constructed using each of the methods allows one to assign an estimated binding energy to a given site.
the correlation of these estimated binding energies with the values of log, kd being the dissociation constant, measured for the seven sites  as described in the materials and methods section is shown in fig.
<dig>  from the summary of the correlations in table  <dig>  it is apparent that the selex procedure produced significantly better information about the binding characteristics of cap than were available using the known sites, while using qpmeme to infer the binding parameters εiα produced better correlation with the measured values of log when used with the selex data set.
the p-values associated with the correlation coefficients  in table  <dig> are as follows: r =  <dig>  => p =  <dig> %, r =  <dig>  => p =  <dig> %, r =  <dig>  => p =  <dig> %.
the best correlation coefficient, obtained for qpmeme trained on selex data, r =  <dig> , represents very significant but less than perfect correlation.
however, one has to remember that the measurements of kd's themselves have a certain amount of error.
therefore, even with perfect predictive power, we would not get a correlation coefficient of one.
note that the number of sequences used from the selex data is comparable to the number of biological cap binding sites used in weight matrix determination.
this observation suggests that the improvement is due to unbiased sampling of binding sequences and not due to the greater number of sites used.
comparison of binding energies of orthologous sites for e. coli and s. typhimurium
without evolutionary pressure to keep the binding energy constant over time, the binding energy of the orthologous site will drift towards the average binding energy, which is set to zero in our convention.
for an e. coli site which is estimated to be a strong binding site, the vast majority of mutations will result in a weaker estimated binding energy.
this has the consequence that, even if the true binding energy is conserved, a poor method of estimating binding energies will probably assign a weaker binding energy to the s. typhimurium orthologs of those e. coli sites with the strongest estimated binding energy.
for each site, s, in e. coli with an s. typhimurium ortholog, s', and for each estimation of the energy matrix, ε, one can define the drift of s according to ε as:
ds=Θe2−e1σe
where e <dig> = ε • s, e <dig> = ε • s' and σe <dig> is the standard deviation of the conditional distribution p .
the function Θ  is the heaviside theta function which  <dig> for x ≥  <dig> and  <dig> for x <  <dig>  the quantity ds measures how much the binding energy of the s. typhimurium site has apparently weakened compared to a measure of how much it would be expected to drift if there were no evolutionary pressure to sustain it.
for each of the four estimations of the energy parameters, the total drift for the highest-scoring  <dig> e. coli binding sites with orthologs in s. typhimurium is shown in table  <dig>  false positives, that is, candidate binding sites identified in e. coli, on the basis of inferred binding energy, which are not functional in vivo, would not be expected to have orthologs in s. typhimurium with significant binding energy.
the overall drift indicated in table  <dig> then receives contributions from both the false positives picked up during the scan of the e. coli genome, and from functional sites whose estimated binding energy differs from the actual binding energy.
figures 2a and 2b show the amounts of drift of individual candidate sites ordered according to drift for the energy parameters inferred from the selex training set and the training set consisting of the known sites, respectively.
overall, the binding parameters inferred using qpmeme systematically indicated less drift than the parameters inferred using the weight matrix, with the best overall performance when the selex data set was used to train the algorithm.
CONCLUSIONS
our purpose in this analysis was to show that indeed with the appropriate kind of training data, one could improve the ability to predict physical and functional binding.
this is in marked contrast to the general feeling in the biology community that many pleoitropic transcription factors bind at too many places and that it is hopeless to try to get functional sites out of motif searches.
recently there were chip-chip experiments done on cap  <cit> .
the conclusion of the study was that cap is physically bound at many thousands of places in the genome.
we saw, however, from the comparative study, that stronger binding sites are significantly conserved, indicating selective pressure.
these results are in agreement with similar studies done in yeast  <cit>  based on dna chip data.
one of the main goals of the experimental procedure was to gather specificity data for a transcription factor  at concentrations comparable to the cellular abundances.
one disadvantage of doing the selection experiment with very high abundance of tf is that it is possible to select sequences where the likely binding site partially overlaps with the primer.
for low abundances the sites bound tend to be in the variable  region, because binding partially to the primer in any window requires enough free energy to be prohibitive.
for high abundance, the threshold for tolerance of such energetic costs is higher.
we believe that, in the original selex sage work  <cit> , the tf abundance is still low enough so that the "primer contamination" problem is avoided.
however, when using this method for a genuinely pleiotropic tf at cellular concentrations, one would face the primer contamination problem we describe here.
one possible resolution of this problem is computational.
we could develop a more complex probability model allowing for the primer contribution and utilize the full data set.
we could do this by generalizing the model in reference  <cit> , and allowing binding in different windows on a longer sequence.
however, to settle the question of whether selex data sets provide any advantage, we decided to focus on small number of sequences where the binding is likely to be in the random or n <dig> part.
the number of selex sequences used for training is comparable to the number of biologically known sites.
hence, the improvement in prediction should be from better sampling of sequences in the dataset rather than from the sample size.
the use of the full data set requires a new computational method, which would be a promising subject for future research.
in this study, we only used sequence data for estimating the parameters related to the motif.
we measured relative affinities  to test the accuracy of our predictions.
as has been shown, combining selex with quantitative affinity measurements leads to even better predictive power  <cit> .
our reason for focusing solely on sequence data is that we foresee developments in short-read sequencing  <cit>  which is expected to lead to a readily available inexpensive technology for generating large selex data sets.
in the last few years, we have seen considerable activity that centers around biophysical aspects of gene regulation.
the push has come from two different ends.
on the one hand, detailed structural modeling of protein-dna interactions has been used to calculate sequence dependent protein dna interaction free energies  <cit> .
the insights from structural considerations could guide the appropriate parametrization of knowledge-based bioinformatic motif discovery tools as well  <cit> .
on the other hand, hidden markov models, which can be thought of as one dimensional statistical mechanics models of multiple proteins binding on dna, have been applied to study multiple binding sites for protein complexes in a stretch of dna  <cit> .
this approach has been extended to nucleosome positioning as well  <cit> .
the remarkable success of these approaches seems to suggest that much could be done with an accurate biophysical description of protein-dna interaction in the context of gene regulation.
precise characterization of the probabilistic protein dna interaction code is a crucial element of such a description.
as technologies for massively parallel signature sequencing  <cit>  become more accurate, large scale selex studies for determining the interaction code would become more and more feasible.

 understanding regulatory circuits controlling gene expression is one of the fundamental problems in modern biology. gene expression is controlled at many different levels but control of transcription is one of the main steps of regulation. one of the best understood control mechanisms is the binding of transcription factors to the regulatory sites on dna in a sequence-specific manner, which affects transcription initiation. the important problem of locating the binding sites for specific tfs, and thus identifying the genes they regulate, has attracted much attention from the bioinformatics community. different methods have been employed for abstracting patterns or "motifs" from the sequences that bind particular tfs leading to predictions of likely binding sites in the genome of the organism under study. factors regulating multiple genes often have binding motifs low in information content, making the task of prediction harder. examples of such highly pleiotropic proteins range from global regulators in prokaryotes to hox proteins, important in metazoan development. experimental approaches to locating binding sites on dna, have uncovered numerous binding sites for various factors. however, looking at the databases devoted to such regulatory sites, like dpinteract and regulondb for e. coli, scpd for yeast and transfac for many higher eukaryotic organisms, it is obvious that, for most pleiotropic tfs targeting a large number of genes, the number of known sites is still a small fraction of all the functional sites. a high-throughput version of the chromatin immunoprecipitation method, commonly known as the "chip on chip", has been introduced recently. in principle, this method locates binding sites genome-wide. however, the resolution is limited to several hundred bases and requires further bioinformatic analysis. an alternative approach would be to find the dna binding specificity of a tf by an in vitro method and then use the binding motif to search the genome for putative sites. one of these methods is selex, which is often used to find the strongest binding sites from a library consisting of randomly generated oligonucleotides. however, a tf can often function at binding sites that are far weaker than the consensus. therefore, to characterize the binding preferences of a tf, we need to identify many of these potential weak binding sites and to estimate the parameters describing the statistical distribution of those sequences. the appropriate modification of the selex procedure needed to achieve this goal is based on the selex-sage procedure. analysis of the conditions under which we get a significant number of intermediate strength sites was performed in. we will use this procedure on the pleiotropic e. coli factor cap. an alternative to this technology would have been to use dna chips for protein binding. currently, for transcription factors with long binding sites, it is common practice to use genomic sequences rather than random libraries in dna chips. this has its advantages but also might lead to uncertainties regarding the genomic background model in the final statistical analysis. to abstract a motif from the sequences found by the modified selex process, we need a computational method: a supervised algorithm, trained on a set of binding sites identified directly by experimental measurements. we will compare different supervised methods for extraction of parameters and use cap targets as a benchmark. the widely used bioinformatic tool for quantitatively describing such motifs is the weight matrix method. setting the threshold correctly is essential for the quality of predictions. however, optimization of the threshold is a non-trivial problem, resolving which is one of the goals of this study. we have shown that using the physically correct expression for binding probability, with saturation effects built in, leads to a more accurate estimate for the binding energy and provides a practically useful solution to the problem of classifier threshold choice. the resulting method, quadratic programming method of energy matrix estimation or qpmeme, turns out to be a one-class support vector machine. in this paper we do the following: perform high throughput, low stringency selex experiments. analyze selex experiments and extract parameters for models of sequence-dependent tf/dna interactions by using qpmeme and its extensions. verify in vitro affinity predicted by models built on selex data by electrophoretic mobility shift assay. identify cap targets in e. coli. compare various predicted sites from different methods to see their functional conservation among gram-negative facultative anaerobes.  sage selex study of cap we followed the selex-sage procedure. after extracting a large number of sequences, we passed these sequences through a quality filter to find sequences that were likely to have been bound by a cap dimer in the random region of the selex pool. we were left with sequences at the end. analysis of sequences from sage selex study the maximum-likelihood method for distributions with sharp cutoffs was described in. in particular circumstances, the parameter estimation method becomes a support vector machine, and the resulting algorithm, qpmeme, can be used to determine the binding energy of the protein to the sequence s, as e = ∑iα εiα siα. the variable siα is defined as follows: if the ith base is α in sequence s, then siα = and siα = otherwise. the parameter εiα is the contribution to the binding free energy from base α being at position i. these εiα parameters are chosen to minimize the variance of e over the background distribution of sequences s, subject to the constraints e) ≤ - for the set of example binding sequences s, j =..., n. sequences satisfying e) ≤ - are then declared to be binding sites. in practice, the base frequencies are taken to be independent and the probability of finding the base α is taken to be pα. the quantity to be minimized is given by ∑iαpαεiα subject to the constraints ∑αpαεiα= for each i. the qpmeme algorithm was used to produce an energy matrix,, using the set of known cap binding sites from the dpinteract database. an energy matrix was also constructed using the binding sites identified by the selex procedure described below. weight matrices were also constructed for both the known sites and the selex sites using the formula wiα = log, where wiα is the i, α component of the weight matrix and fiα=∑jsiα is the frequency of the ith base being α. the background probabilities for g and c are taken to be the same. the same applies to the background probabilities for a and t. for both algorithms, the background gc content was taken to be when constructing energy matrices from the known sites, and was taken to be when constructing energy matrices from selex sites. for the weight matrices/energy matrices used in the study, see additional file additional file additional file and additional file prediction of in vitro affinity of binding sites the energy matrices constructed using each of the methods allows one to assign an estimated binding energy to a given site. the correlation of these estimated binding energies with the values of log, kd being the dissociation constant, measured for the seven sites as described in the materials and methods section is shown in fig. from the summary of the correlations in table it is apparent that the selex procedure produced significantly better information about the binding characteristics of cap than were available using the known sites, while using qpmeme to infer the binding parameters εiα produced better correlation with the measured values of log when used with the selex data set. the p-values associated with the correlation coefficients in table are as follows: r = => p = %, r = => p = %, r = => p = %. the best correlation coefficient, obtained for qpmeme trained on selex data, r =, represents very significant but less than perfect correlation. however, one has to remember that the measurements of kd's themselves have a certain amount of error. therefore, even with perfect predictive power, we would not get a correlation coefficient of one. note that the number of sequences used from the selex data is comparable to the number of biological cap binding sites used in weight matrix determination. this observation suggests that the improvement is due to unbiased sampling of binding sequences and not due to the greater number of sites used. comparison of binding energies of orthologous sites for e. coli and s. typhimurium without evolutionary pressure to keep the binding energy constant over time, the binding energy of the orthologous site will drift towards the average binding energy, which is set to zero in our convention. for an e. coli site which is estimated to be a strong binding site, the vast majority of mutations will result in a weaker estimated binding energy. this has the consequence that, even if the true binding energy is conserved, a poor method of estimating binding energies will probably assign a weaker binding energy to the s. typhimurium orthologs of those e. coli sites with the strongest estimated binding energy. for each site, s, in e. coli with an s. typhimurium ortholog, s', and for each estimation of the energy matrix, ε, one can define the drift of s according to ε as: ds=Θe2−e1σe where e = ε • s, e = ε • s' and σe is the standard deviation of the conditional distribution p. the function Θ is the heaviside theta function which for x ≥ and for x < the quantity ds measures how much the binding energy of the s. typhimurium site has apparently weakened compared to a measure of how much it would be expected to drift if there were no evolutionary pressure to sustain it. for each of the four estimations of the energy parameters, the total drift for the highest-scoring e. coli binding sites with orthologs in s. typhimurium is shown in table false positives, that is, candidate binding sites identified in e. coli, on the basis of inferred binding energy, which are not functional in vivo, would not be expected to have orthologs in s. typhimurium with significant binding energy. the overall drift indicated in table then receives contributions from both the false positives picked up during the scan of the e. coli genome, and from functional sites whose estimated binding energy differs from the actual binding energy. figures 2a and 2b show the amounts of drift of individual candidate sites ordered according to drift for the energy parameters inferred from the selex training set and the training set consisting of the known sites, respectively. overall, the binding parameters inferred using qpmeme systematically indicated less drift than the parameters inferred using the weight matrix, with the best overall performance when the selex data set was used to train the algorithm.  our purpose in this analysis was to show that indeed with the appropriate kind of training data, one could improve the ability to predict physical and functional binding. this is in marked contrast to the general feeling in the biology community that many pleoitropic transcription factors bind at too many places and that it is hopeless to try to get functional sites out of motif searches. recently there were chip-chip experiments done on cap. the conclusion of the study was that cap is physically bound at many thousands of places in the genome. we saw, however, from the comparative study, that stronger binding sites are significantly conserved, indicating selective pressure. these results are in agreement with similar studies done in yeast based on dna chip data. one of the main goals of the experimental procedure was to gather specificity data for a transcription factor at concentrations comparable to the cellular abundances. one disadvantage of doing the selection experiment with very high abundance of tf is that it is possible to select sequences where the likely binding site partially overlaps with the primer. for low abundances the sites bound tend to be in the variable region, because binding partially to the primer in any window requires enough free energy to be prohibitive. for high abundance, the threshold for tolerance of such energetic costs is higher. we believe that, in the original selex sage work, the tf abundance is still low enough so that the "primer contamination" problem is avoided. however, when using this method for a genuinely pleiotropic tf at cellular concentrations, one would face the primer contamination problem we describe here. one possible resolution of this problem is computational. we could develop a more complex probability model allowing for the primer contribution and utilize the full data set. we could do this by generalizing the model in reference, and allowing binding in different windows on a longer sequence. however, to settle the question of whether selex data sets provide any advantage, we decided to focus on small number of sequences where the binding is likely to be in the random or n part. the number of selex sequences used for training is comparable to the number of biologically known sites. hence, the improvement in prediction should be from better sampling of sequences in the dataset rather than from the sample size. the use of the full data set requires a new computational method, which would be a promising subject for future research. in this study, we only used sequence data for estimating the parameters related to the motif. we measured relative affinities to test the accuracy of our predictions. as has been shown, combining selex with quantitative affinity measurements leads to even better predictive power. our reason for focusing solely on sequence data is that we foresee developments in short-read sequencing which is expected to lead to a readily available inexpensive technology for generating large selex data sets. in the last few years, we have seen considerable activity that centers around biophysical aspects of gene regulation. the push has come from two different ends. on the one hand, detailed structural modeling of protein-dna interactions has been used to calculate sequence dependent protein dna interaction free energies. the insights from structural considerations could guide the appropriate parametrization of knowledge-based bioinformatic motif discovery tools as well. on the other hand, hidden markov models, which can be thought of as one dimensional statistical mechanics models of multiple proteins binding on dna, have been applied to study multiple binding sites for protein complexes in a stretch of dna. this approach has been extended to nucleosome positioning as well. the remarkable success of these approaches seems to suggest that much could be done with an accurate biophysical description of protein-dna interaction in the context of gene regulation. precise characterization of the probabilistic protein dna interaction code is a crucial element of such a description. as technologies for massively parallel signature sequencing become more accurate, large scale selex studies for determining the interaction code would become more and more feasible. 
the consequences of this increase in accuracy for prediction of in vivo binding sites in the e. coli genome are also discussed. characterizing transcription factor binding motifs is a common bioinformatics task. for transcription factors with variable binding sites, we need to get many suboptimal binding sites in our training dataset to get accurate estimates of free energy penalties for deviating from the consensus dna sequence. we believe that with progress in short-read sequencing technology, one could use selex methods to characterize binding affinities of many low specificity transcription factors. to obtain large number of sequences required for this analysis we used a selex sage protocol developed by roulet et al. we found that bioinformatics scores based on selex sage data does better in terms of prediction of physical binding energies as well as in detecting functional sites. we analyzed low stringency selex data for e. coli catabolic activator protein, and we show here that appropriate quantitative analysis improves our ability to predict in vitro affinity.
 characterizing transcription factor binding motifs is a common bioinformatics task. for transcription factors with variable binding sites, we need to get many suboptimal binding sites in our training dataset to get accurate estimates of free energy penalties for deviating from the consensus dna sequence. one procedure to do that involves a modified selex method designed to produce many such sequences.  we analyzed low stringency selex data for e. coli catabolic activator protein, and we show here that appropriate quantitative analysis improves our ability to predict in vitro affinity. to obtain large number of sequences required for this analysis we used a selex sage protocol developed by roulet et al. the sequences obtained from here were subjected to bioinformatic analysis. the resulting bioinformatic model characterizes the sequence specificity of the protein more accurately than those sequence specificities predicted from previous analysis just by using a few known binding sites available in the literature. the consequences of this increase in accuracy for prediction of in vivo binding sites in the e. coli genome are also discussed. we measured the dissociation constants of several putative cap binding sites by emsa and compared the affinities to the bioinformatics scores provided by methods like the weight matrix method and qpmeme trained on known binding sites as well as on the new sites from selex sage data. we also checked predicted genome sites for conservation in the related species s. typhimurium. we found that bioinformatics scores based on selex sage data does better in terms of prediction of physical binding energies as well as in detecting functional sites.  we think that training binding site detection algorithms on datasets from binding assays lead to better prediction. the improvements in accuracy came from the unbiased nature of the selex dataset rather than from the number of sites available. we believe that with progress in short-read sequencing technology, one could use selex methods to characterize binding affinities of many low specificity transcription factors. 