BACKGROUND
the specific recognition of genomic cis-regulatory elements by nucleic acid binding proteins is of critical importance for many vital processes such as dna replication and repair, mrna translation and transcriptional regulation.
specific tf-dna interactions, as well as protein-protein contacts established between closely located tfs at promoter regions, are key determinants of the concerted expression of genes in response to external stimuli.
the task of determining the relationships among tfs and their targets has been widely addressed experimentally  <cit>  and computationally  <cit> .
of great relevance and use in recent years are computational approaches, based on the construction of statistical models that characterize the dna-binding preferences of tfs, which are in turn used to scan genomic sequences in an effort to identify new putative binding sites  <cit> .
the probabilistic models more commonly used in such computational approaches are position weight matrices  obtained from multiple alignments of known binding sites.
this approach is limited to tfs with a sufficient number of experimentally identified binding sites, for which reliable statistical models may be built.
an alternative approach would be to predict dna operator sites which are compatible with the mode of binding of a given tf  <cit> .
this approximation would depend on two key components: i) the knowledge of the protein and dna residue positions involved in binding at the spatial level and ii) a method to evaluate the compatibility of different dna bases and amino acids to interact  <cit> .
the computational analysis of available protein-dna complexes has resulted in several papers describing the characteristics of the amino acid-base interactions that determine binding specificity  <cit> , the different types of readout mechanisms involved in dna recognition  <cit>  and the evolutionary conservation of residues located at contact interfaces  <cit> .
direct readout is associated to recognition through contacts established between atoms from amino acid side-chains and nitrogen bases.
indirect readout, on the other hand, is mediated by the contribution of residues of the proteins and dna which are not in direct contact and conformational changes undergone by dna upon protein binding  <cit> .
these reports show that both direct and indirect readouts significantly contribute to specific protein-dna recognition.
with respect to direct readout, it is estimated that about two-thirds of all protein-dna interactions are van der waals contacts which do not generally confer sequence specificity, with the exception of the hydrophobic interactions involving the c <dig> atom of thymine  <cit> .
on the contrary, hydrogen bonds are the major source of specific interaction contacts.
two-thirds of the hydrogen bonds formed between amino acids and bases are bi-dentate or complex interactions providing a great specificity.
non-classical c-h ···o hydrogen mediated links have also been found at protein-dna contact interfaces  <cit> , albeit their energetic contributions to recognition are not fully understood yet.
these interactions are inherently weak individually  <cit> , which means that cumulative effects are indispensable to significantly account for specificity, as occurs with hydrophobic links.
water mediated bridges, though common in interaction interfaces, are mostly used as gap-fillers  <cit>  and are also engaged in buffering unfavorable electrostatic repulsions between interacting atoms at the interface  <cit> .
electrostatic interactions, most of which are established between the protein main chain and the sugar/phosphate backbone, do not significantly contribute to specific recognition, although they play important roles in the transition from the unspecific to the specific complex  <cit> .
to date just a few reports have been published aiming at finding putative tf target sequences using structural information.
there have been some attempts to apply physical potential functions – i.e., in the form of atomic force fields – to estimate the energy of interaction and the relative contribution of direct and indirect readout mechanisms  <cit> .
however, these physical models do not appear to significantly outperform simpler statistical methods  <cit> .
examples of these methods are approaches oriented to extracting structural information from family-wise comparisons, building a statistical model of the interface derived from crystallographic structures and binding sites of proteins that belong to the same family  <cit> .
other procedures infer statistical potentials of interaction at the residue-base level from datasets of known complexes  <cit> , which are used to examine the compatibility between the protein and its putative sites.
however, there is some evidence suggesting that the performance of these methods may be limited by the simplicity of their interaction potentials, calculated at the level of cα atoms for protein residues.
moreover, the increasing number of high-resolution structures for protein-dna complexes opens a new door towards a more detailed study of their contact interfaces.
indeed a recent report confirms that atomic details improve the ability of physical potential functions to predict sequence-specific protein-dna interactions  <cit> .
in this work we constructed position weight matrices that capture the binding specificity of transcription factors, based on information extracted from the protein data bank .
three atomic preference matrices, for hydrogen bonds, water-mediated hydrogen bonds and hydrophobic interactions, were derived from a non-redundant training set of  <dig> complexes annotated in the pdb  <cit> .
these matrices were used to make explicit atomistic representations of hydrogen and hydrophobic bonds as well as the contribution of water molecules at interfaces, which gives us the opportunity to score the direct readout.
this contribution is combined with empirical estimations of dna deformation, in order to calculate a potential of binding that includes both direct and indirect readout contributions.
we evaluated the performance of our algorithm in a set of  <dig> bacterial and  <dig> eukaryotic tfs which have been co-crystallized bound to dna, and in most cases our results proved to be as good as or better than those obtained with the structure-based cumulative contact method by morozov and siggia  <cit> .
in addition, two tf homology models were analyzed in detail and used to predict their dna binding motif, after sampling side-chain rotamers at their contact interfaces.
in this case the results we obtained were significantly better than those returned by the reference method, which indicates that our algorithm could be suitably used to study tfs of unknown structure starting from structural models.
we also discuss the strengths and limitations of our approach that might potentially be used for tfs with few or no experimentally characterized binding sites.
RESULTS
protein-dna interface atomic contacts and interaction preferences
starting from a non-redundant training set of crystallographic protein-dna complexes culled from the pdb, we constructed atomic frequency matrices for hydrogen bonds, water-mediated hydrogen bonds and hydrophobic contacts.
a close inspection at the information embodied in these frequency matrices  reveals some interesting features of protein-dna interfaces.
for instance, as claimed in previous reports, we found that arginine is the major source of hydrogen bonds, with a marked preference to interact with guanine  <cit> .
however, in the present study we divided this preference in pairs of interacting atoms, finding that groups nh <dig> and nh <dig> of arginine establish 88% of the hydrogen bonds with atoms n <dig> and o <dig> from guanine.
histidine also showed a marked preference towards guanine, with atom ne <dig> interacting in  <dig> out of  <dig> hydrogen bonds with atoms o <dig> and n <dig> from the nitrogen base.
overall, we found  <dig> interface hydrogen bonds in our training set.
hydrogen bonds reported by hbplus  <cit>  between pairs of atoms in the side-chains of amino acids and bases in the set of  <dig> protein-dna complexes extracted from pdb were computed.
rows correspond to the hydrogen bond donor/acceptor atoms from the amino acids and columns to the hydrogen bond donor/acceptor atoms from bases  with their corresponding pdb atom codes.
each cell indicates the number of contacts found for a given pair of atoms.
in contrast with hydrogen bonds, the landscape of hydrophobic interactions is quite different.
as previously reported  <cit> , the main discriminatory group regarding specific recognition through hydrophobic interactions is the methyl group of thymine.
accordingly, the c <dig> group accounted for 20% of all  <dig> hydrophobic interactions found in our training dataset, being the main source of contacts for all amino acids when compared with the other nitrogen bases.
with the exception of c <dig>  we found no obvious interaction preferences between amino acid and bases.
this means that the signal to noise ratio of the calculated preferences is lower than the one observed for hydrogen bonds, and are therefore less informative.
thus, in our subsequent analyses, we decided to take into account only the contribution of the methyl group of thymine during the evaluation of the interface between a tf and a given dna sequence.
we also generated frequency matrices for water-mediated hydrogen bonds, finding a total of  <dig> atomic interactions.
the overall observed interaction propensities are quite similar to those in the hydrogen bond matrix, which is probably related to the gap filling role of water that supports the existence of long-distance hydrogen bonds  <cit> .
in addition, the buffering role of water allows the formation of electrostatically unfavorable hydrogen bonds, adding new interaction propensities absent in the hydrogen bond matrix.
for instance, we found eight water-mediated hydrogen bonds between arginine side-chain nitrogen groups and cytosine n <dig>
a survey of the quality of the atomic interaction matrices
the redundancy of the data in the pdb and the bias towards specific protein folds considerably tangle the efforts aiming at exploiting the wealth of this database, limiting the scope of the results.
in order to overcome these problems, we decided to avoid redundancy as much as possible in our training set, taking care of withholding a sufficiently informative set of pdb entries.
we also planned a thorough bootstrap assay to estimate the quality of the atomic preferences extracted from the training dataset
in this bootstrap tryout we resampled with replacement a thousand subsets of entries, randomly excepting a considerable fraction of the entries included in the initial training set as described in methods.
by leaving out almost a third of the total entries, we measured the degree of overtraining due to residual redundancy in the initial set, as well as the statistical relevance of the interaction preferences reckoned in the matrix building process.
a straightforward way of addressing this issue is by finding out whether the matrices constructed using the truncated training sets were good enough to suitably evaluate the direct readout for entries used to build the corresponding bootstrap matrices, though not failing to score the crystallographic complexes that were excluded in the resampling step.
as may be seen in figure  <dig>  the exclusion of 30% of the entries in the training set does not significantly affect the discrimination capability of the matrices.
here we assess the scores obtained when evaluating the bootstrapped  and excluded  datasets with the bootstrap matrices – in the abscissa -, against the generic  or shuffled matrices .
we found in the scatter plots and the regression analysis that the scores obtained for the assay of the bootstrapped and excluded datasets with the bootstrap matrices correlate very well – i.e., r <dig> =  <dig>  and r <dig> =  <dig>  respectively -, with the values obtained when assaying those same datasets with the generic matrices .
as anticipated, the correlation entirely disappears when the analysis is done with the shuffled matrices – two bottom charts of figure  <dig>
direct and indirect readout mechanisms in protein-dna complexes
the information embodied in the atomic interaction matrices presented in the previous section can be used to estimate the contribution of direct readout in protein-dna recognition.
however, it has been shown that a mechanism of indirect readout also plays a relevant role in this process, at least for some tfs, and was therefore considered in this work.
in particular, we modeled indirect readout as the cost of threading a nucleotide sequence into a fixed dna backbone.
in order to integrate both direct and indirect recognition mechanisms, we designed a saturating mutation strategy, which is the kernel of the dnaprot algorithm.
briefly, the algorithm iteratively evaluates the interaction potential of a given tf as the docked nucleotide sequence mutates in a 4n space, for a dna duplex of length n .
this renders a structure-based position weight matrix that can be used to scan genomic sequences.
a full description is provided in the methods section, but it is important to note that direct and indirect readout scores are linearly combined by means of a deformation weight, d. as d gets larger, the relevance of indirect readout increases.
a study of the performance of our method is shown in figure  <dig>  where it can be observed that the indirect readout contribution to the ability of correctly scoring cognate sites of crp and narl is critical, since the predictive competence increases with large d values, as revealed by the increasing areas under the curve  in the roc plots.
these outcomes agreed with previous reports in the literature claiming the central role of deformation energy in the site recognition mechanism of crp  <cit>  and narl  <cit> , which bend dna 90° and 42° respectively.
purr and dnaa display an opposite trend, since their largest aucs are obtained with d values less than  <dig> .
while purr is known to bend dna with an angle of 45°  <cit> , the reported bending angle of dnaa is just 28°  <cit> .
this means that although in the first two cases the contribution of deformation can be considered essential, in the cases of purr and dnaa optimal predictive values are obtained with intermediate values of d, also including an important contribution of direct readout.
the subtle changes observed for the auc in the roc plots of figure  <dig>  in fact correspond to significant variations in the classification performance of our method.
this can be seen in the results reported in table  <dig> for the scoring of purr known binding sites.
the comparison of the pwms obtained with d values of  <dig>  and  <dig> , those for which we obtained the best and worst performance in figure 3c, reveals that setting d to its optimum value causes the p-values of cognate sites to be reduced by one or two orders of magnitude.
this improvement corresponds to a reduction of 103- <dig> false positives extracted along with true sites in a search set of 106bp.
it is worth noting that, using the optimum d value in the construction of the model for this tf and setting the sensitivity to 80%, the number of false positive sites recovered with the true sites is as low as  <dig> in a genome-sized search set.
the known binding sites of purr reported in regulondb  <cit>  were scored using dnaprot, maintaining w fixed to  <dig> and setting d to  <dig>  and  <dig> .
these values of d correspond to the parameterization that renders the worst and the best pwms for this transcription factor respectively.
the first column corresponds to the binding sites, the second and fourth columns report the scoring obtained with dnaprot and the third and fifth columns show the p-values of the scores obtained in a site recovery assay in a background set of  <dig> random sequences.
scoring crystallographic protein-dna complexes with the dnaprot algorithm
in this section we present the results of scoring a diverse collection of protein-dna complexes solved by x-ray crystallography.
the test set includes eight transcription factors bound to cognate dna sequences, of which four are prokaryotic and the other four eukaryotic, belonging to eight different superfamilies according to scop  <cit> .
in this experiment we took the coordinates of each of these interfaces and compute structural position weight matrices  that approximate their binding specificity.
two types of pwms are calculated here: i) readout matrices, derived using the dnaprot algorithm outlined in figure  <dig>  and ii) cumulative contact matrices, obtained by simply counting the contacts between protein side-chains and nitrogen bases, as previously described by morozov and siggia  <cit> .
these pwms can then be compared with cognate matrices built from experimentally determined binding sites, by means of local alignments  <cit> .
further details are provided in the methods section, but it is relevant to note that readout pwms are computed using knowledge-based atomic potentials, whereas cumulative contact pwms are calculated with the assumption that the consensus dna sequence is part of the pdb complex.
figure  <dig> includes the sequence logos derived from the computed pwms; the expectation values of pwm comparisons are shown in table  <dig>  together with details of the test set.
scop v <dig>
in the first column we include the name of the structural superfamily of each tf according to scop  <cit> .
in the second column the name of each tf is followed by the pdb identifier of the structure used enclosed in brackets.
the following two columns contain the resolution and the r-value of the crystals respectively.
the last two columns report the e-values resulting from stamp  <cit>  local alignments between the structure-based pwms obtained with the method of morozov and siggia  <cit>   or our approach  and the corresponding cognate pwms, extracted from the literature  <cit> .
parenthesized e-values were calculated after sampling interface rotamers with dnaprot.
the pwms generated by dnaprot starting from crystallographic structures are comparable or better than those produced with the methodology of cumulative contacts by morozov and siggia  <cit> .
this is clear in the results of table  <dig>  where, with the exception of purr and leu <dig>  the statistical significance of our readout matrices is higher or of the same order as that of the contact matrices.
an inspection of figure  <dig> also shows that for the three prokaryotic regulatory proteins, crp, narl and metj, the sequence logos obtained by both methods are quite similar.
however, except for the case of metj, for which the significance of the pwms generated by both methods is of the same order of magnitude, for crp and narl the matrices generated by our algorithm are two orders of magnitude better.
on the contrary, although the logo of purr is correct, the e-value it gets when compared with the expected logo is several orders of magnitude worse, as a consequence of failing to identify  <dig> consensus nucleotides.
the results for our selection of eukaryotic tfs show close similarities for the statistical significance of the pwms generated by both methods.
however, for rap <dig> sampling of interface rotamers caused an improvement of the readout matrix that outperforms the contact matrix by one order of magnitude, .
something similar happened for the matrix generated for purr, though in this case the improvement of three orders of magnitude of the significance is still worse than the value of the contacts pwm.
footprinting comparative models of protein-dna complexes with dnaprot
in this section we further test the ability of dnaprot readout matrices to capture the binding specificity of tf-structural models.
here we generated homology models for two transcription factors: fnr, a global regulatory protein in e. coli, and giant, a regulatory protein involved in the early development of d. melanogaster.
as in the previous section, we also computed cumulative contact pwms for comparison purposes, and we evaluated both structure-based matrices by aligning them to cognate matrices built from operator sites reported in the literature.
we now present the results of both models in more detail.
the protein sequence of fnr can be confidently aligned to the sequence of e. coli crp  with tfmodeller  <cit> , despite an overall sequence identity of just 22%, providing a reliable frame for comparative modeling.
however, several residue substitutions occur at the binding interface, suggesting that the readout might be changing from the template to the model.
indeed we found that four out of ten interface residues are mutated in the dimeric model, as annotated in table  <dig>  for this reason we sampled different rotamers for the side-chains of the interface residues and selected those that yielded best readout scores using the dnaprot algorithm.
as shown in figure  <dig>  this limited sampling was enough to improve the readout pwm, increasing the number of correct consensus nucleotides in the sequence logo from six to nine.
this improvement of the model is also related to an increase of the statistical significance of the readout matrix of one order of magnitude with respect to the cumulative contact matrix, as shown in table  <dig>  a further improvement is also possible by deriving a contact matrix from the refined homology model, which causes an additional increase of the e-value to  <dig> e- <dig>  the multiple alignment, figure 5d, also suggests that the interface differences inferred from the comparative model correlate with the sequence conservation of the dna-binding domain among fnr and crp homologous sequences.
scop v <dig>
winged helix
leucine zipper domain
in the first column we include the name of the structural superfamily of each tf according to scop  <cit> .
the name of each model is followed by the pdb identifier of the template used to guide model building, enclosed in brackets.
the following three columns contain the ratio of contact interface identity, the resolution and the r-value of the crystals respectively.
the last two columns report the e-values resulting from stamp  <cit>  local alignments between the structure-based pwms obtained with the method of morozov and siggia  or our approach  and the corresponding cognate pwms, extracted from the literature  <cit> .
parenthesized e-values were calculated after sampling interface rotamers while bracketed e-values correspond to cumulative contact pwms computed over the refined dnaprot complex.
the second study of a homology model was performed with giant, an eukaryotic leucine zipper transcription factor that has been studied previously  <cit> .
this protein can be modeled based on human cebpb , with a sequence alignment that covers 92% of the dna-binding domain.
although the overall sequence identity is only 13%, six out of  <dig> interface residues, as labeled by tfmodeller, are conserved.
the first evaluated giant model yielded a pwm that captures only four correct consensus positions, as seen in figure  <dig>  as in the case of fnr, we refined the binding interface of giant by sampling rotamers and selecting those with best readout scores.
the refined model yields a structure-based pwm that includes eight consensus nucleotides, improving the preliminary model results in two orders of magnitude in terms of e-value .
this sampled model also suggests a different set of interface residues responsible for dna sequence discrimination.
the figure shows an alignment of giant and cebpb homologous sequences that seem to cluster in two different modes of binding to dna.
relative contribution of atomic interactions at dna recognition interfaces
the dnaprot algorithm models the direct readout mechanism by labeling hydrogen bonds, hydrophobic interactions and water-mediated hydrogen bonds at the interface, which are evaluated by means of log-likelihood preference matrices.
therefore, the data generated for the eight tfs presented in table  <dig> can be further analyzed with the aim of calculating the relative contribution of each of these atomic interactions to dna recognition.
as shown in table  <dig>  we found four cases in which recognition is mediated exclusively by hydrogen bonds , while the remaining tfs seem to employ two or more types of atomic interactions to drive dna binding.
of these, perhaps the most interesting cases are narl, which has a large hydrophobic contribution, and zif <dig>  with 30% of the interface score contributed by water-mediated hydrogen bonds.
repeating these calculations on the set of non-redundant complexes, we found that 78% of dnaprot interface scores are contributed by hydrogen bonds, while water-mediated hydrogen bonds account for 16% and hydrophobic interactions are responsible for 6% of the total.
table 5b shows a more detailed analysis of narl and zif <dig>  in which structure-based pwms derived from only one type of atomic interactions are compared with the same cognate pwms in table  <dig>  these last results show that the combination of explicit representations of the relative contribution of the different types of interaction improves the predictive competence of dnaprot.
a) transcription factors in table  <dig> where further analyzed by assessing the fraction of the interface binding score, obtained after summing all log-likelihood atomic scores, contributed by each interaction type.
the number of atomic interactions found by the dnaprot protocol in the native pdb structure is indicated in the rightmost column.
b) a further analysis of two notable examples, narl with an important hydrophobic contribution and zif <dig> with a remarkable contribution of water-mediated hydrogen bonds, in which structurally derived pwms that consider exclusively h-bonds, hydrophobic contacts or water-mediated h-bonds are compared with default dnaprot pwms, which integrate all three atomic interactions.
discussion
the field of computational prediction of regulatory sites in genomic sequences has been dominated by methodologies relying on prior knowledge of experimentally characterized binding sites  <cit> .
only a limited group of studies has intended to generate structural representations of the sequence recognition interface, using approaches that somehow approximate the physical processes involved in dna recognition  <cit> .
some of these studies have dissected binding interfaces at the residue level  <cit> , while others use atomic-based descriptions  <cit> .
though some efforts employ force fields for the calculations, most of these methodologies rely on pre-compiled interaction preferences.
moreover, there are also conceptually simpler approaches that rely on counting contacts at the interface  <cit>  that generally work under the assumption that the nucleotide sequence captured in a crystallographic protein-dna complex is the consensus.
these contact-based approaches, particularly the cumulative-contact method of morozov and siggia, promise to be successful in binding site searches at the genomic scale if the structure of a bound tf is at hand.
extraction of protein-dna atomic interaction preferences from the pdb
the dnaprot algorithm introduced here incorporates an atom-based representation of protein-dna interfaces and explicitly integrates both direct and indirect readout.
though the thermodynamic contribution of interface desolvation and conformational entropy of interacting groups have been considered in some approaches using complex formulations  <cit> , our algorithm, outlined in figure  <dig>  is concerned with binding specificity and does not account for those contributions, like other related reports  <cit> .
however, our methodology attempts to approximate these phenomena by including explicit representations of binding interface water molecules and side-chain rotamer sampling of interacting amino acids, which have important implications in our results.
the hydrogen bond and hydrophobic interaction matrices used by dnaprot were derived from a non-redundant set of  <dig> complexes, naturally surpassing the training sets used in preceding studies, as more structural data is now available at the pdb.
an important assumption behind the construction of this data set is that rules governing specific recognition of dna by proteins are generally due to the conformational restrictions imposed by the double helix.
therefore, we chose to collect a comprehensive training set instead of using a superfamily-specific or a tf-exclusive set, in which the scantiness of information would have weakened the derived statistical models.
in addition, we built an atomic contact matrix that explicitly accounts for water-mediated hydrogen bonds in protein-dna interfaces, which constitutes a novel contribution of this work.
the data obtained in the construction of our contact propensity matrices, such as the hydrogen bond matrix in table  <dig>  constitute an update of the pioneer work by mandel-gutfreund et al.
<cit>  and luscombe et al.
<cit>  as we found similar trends in the atomic interaction propensities in our enlarged dataset.
in addition, this work presents quantitative arguments , which demonstrate the importance of including hydrophobic and water-mediated hydrogen bonds in the list of interactions that contribute to specific recognition of nucleotide sequences.
assessment of the quality of the atomic interaction matrices
keeping in mind the problem of data bias in the pdb, we planned an in-depth bootstrap tryout to assess the reliability of the information contained in our contact matrices.
the exclusion of almost a third of the entries in the original training set yielded atomic interaction matrices with only a slightly reduced ability to correctly evaluate infrequent atomic contacts.
these results confirm that the preference matrices derived in this work are biologically significant and that the atomic associations found are not significantly affected by overtraining nor by redundancy, since bootstrap matrices scored the excluded training sets almost as well as the bootstrapped ones.
on the contrary, the results for the negative control obtained while repeating this analysis with shuffled matrices , show completely scattered distributions with no correlation.
these outcomes support the consistency of the interaction preferences extracted from the training set used to build our matrices.
they also corroborate that the inclusion of new protein-dna complexes, regularly deposited in the pdb, will not cause considerable variations of the interaction preferences, except for very unusual atomic contacts.
evaluation of direct/indirect readout contributions to protein-dna recognition
the rationale for using atomic potentials is that many interface contacts are either complex or bidentate  <cit> , and therefore cannot be properly accounted for by residue-based approaches.
in those cases the contact pairing ratio is not 1: <dig>  as the same amino acid might be in contact with two base steps simultaneously or binding a given nitrogen base through multiple groups.
these drawbacks could be resolved by using explicit atomistic representations of the mode of binding, which have been the goal of a group of papers in this field  <cit> .
however, using atomic potentials requires high-quality atomic structures of protein complexes.
in addition, atomic-detail approaches such as dnaprot might be more affected by cases in which side-chains rearrange upon mutation of the bound nucleotides  <cit> , which is the reason that it might be necessary to sample interface rotamers while doing comparative modeling, as explained in the results section.
further sampling, including backbone movements, might be required in certain cases, as already envisaged by havranek et al.
<cit> .
in fact, previous works have already recognized that homologous protein-dna interfaces change their docking geometry as their sequences diverge  <cit> .
besides the contribution of the direct readout to sequence recognition, we also considered in our model a mechanism of indirect readout, as sequence-specific dna deformation has been identified as key to dna recognition for many transcription factors  <cit> .
the algorithm presented in this work follows previous efforts that model indirect readout as the cost of deforming a dna duplex  <cit> .
both readout mechanisms, shown to be critical for specific recognition in figure  <dig>  are linearly integrated into a single binding potential.
the weights of both direct and indirect readout can be tuned for different transcription factors according to their docking mode, as the examples depicted in figure  <dig> imply.
this observation suggests that the performance of the dnaprot protocol can be improved using previous biological knowledge and this is certainly a desirable property.
this data also insinuates that each tf might have its own balance between direct and indirect readout, although we cannot exclude the possibility that this value is a property of its structural superfamily.
moreover, data in table  <dig> gives further insight into the impact that correct weighting could have on genomic scale tf-binding-site prediction assays, as the number of false positives might be considerably reduced in order to obtain a more reliable set of predictions.
assaying the predictive potential of dnaprot in crystallographic structures and homology models
the predictive power of dnaprot, evaluated in tables  <dig> and  <dig> and figures  <dig> and  <dig>  suggests that our readout pwms are, with one exception, comparable to or better than those generated by a well established reference methodology, the cumulative contact pwm proposed by morozov and siggia  <cit>  and related to previous reports  <cit> .
it is worth remembering the reader that, unlike the reference method, dnaprot does not assume that the nucleotide sequence in the input pdb complex is the consensus; rather, it performs an in silico mutagenesis assay and evaluates 4n sequences using a readout function.
as already mentioned, the conceptually simpler cumulative contact approach does not consider the contribution of indirect readout.
this might explain the better results obtained with our method for crp and narl shown in table  <dig>  two tfs known to have an important contribution of indirect readout  <cit> .
in contrast, the example of purr shows that the dnaprot formulation of hydrogen bonds and hydrophobic interactions does not fully capture the array of atomic contacts at the interface, as the obtained readout sequence logo misses  <dig> consensus positions that might correspond to a different type of interaction.
despite this fact, the predicted pwm correctly identifies the remaining consensus positions except one, and the statistical significances obtained in the site recovery assay of the cognate binding sites of this tf depicted in table  <dig> prove the predictive competence of the matrix built.
notwithstanding the promising results obtained for crystallographic structures, the most valuable application of the methodology presented in this paper is found in the exercise of comparative modeling.
two examples were modeled here, a well-characterized prokaryotic regulator  and an eukaryotic factor previously studied in a related article  <cit> .
both examples demonstrate that rotamer sampling at the interface is necessary for obtaining optimal results, and we found that the readout formulation presented here can be effective in selecting the best rotamers.
as shown in table  <dig> and figure  <dig>  our methodology yielded better results than those obtained with the cumulative contact method for two homology models.
the reason for this could be that comparative models usually contain errors in the assigned position of protein side-chains and, most importantly, they do not necessarily contain the consensus dna sequence.
a combination of in silico saturating dna mutagenesis and interface side-chain sampling allows the dnaprot algorithm to partially overcome these problems.
not only is sampling positive in modeling tasks, but we also found that the crystallographic structures of purr and rap <dig> yielded better readout pwms after resampling their interface side-chain rotamers.
although the quality of the structures used by our methodology is of primary relevance, we could not find a clear correlation between r-factor, resolution or experimental techniques for entries in the pdb and the outcomes of our procedure.
nevertheless, we often encountered protein-dna complexes where only a few hydrogen bonds or hydrophobic interactions can be identified using standard bond geometries, yielding only partial dna motifs.
in such cases we found that the cumulative contact method seems to be less sensitive to the quality of structural data.
further work needs to be done to explore whether side-chain sampling, and even backbone sampling, can help circumvent the limitations that data quality imposes on the performance of readout pwms.
evaluation of the contribution of different interaction types to the recognition process
the atomistic foundation of our methodology also gives us the possibility of exploring the relative contribution of each interaction type to the dna recognition process.
the analysis of a non-redundant set of more than  <dig> complexes proves that water-mediated hydrogen bonds are the second source of specific interactions at contact interfaces, which raises a warning regarding the usual exclusion of these interactions in structural studies trying to model protein-dna interaction interfaces.
the exclusion of water overlooks a wide group of highly informative contacts, mainly long-distance and unfavorable hydrogen bonds  <cit>  that constitute novel pairings absent from the hydrogen bond preference matrix.
the example of zif <dig> analyzed in table 5b, shows that an important part of the information content of the pwm obtained for this tf corresponds to water-mediated interactions, and the inclusion of this information considerably increases the statistical significance of the readout matrix.
something similar occurs for hydrophobic interactions; despite being relatively infrequent sometimes their contribution to binding might be central.
this is the case of the matrix obtained for narl, in which hydrophobic links account for most of the information content of the structurally-derived pwm, despite being the least frequent interface atomic interaction.
in this last example, without a reliable inclusion of those infrequent interactions, our model for this tf would have been worthless.
CONCLUSIONS
in summary, our results suggest that the dnaprot algorithm, together with the set of atomic interaction matrices obtained in this work, have useful applications.
the matrices contain biologically meaningful information that confirm and enrich previous reports at the atomic level of interaction.
in addition to the uses presented in this paper, these matrices could be taken to estimate specificity of binding  <cit>  or as a guide for crystallographic studies of protein-dna complexes.
with respect to the algorithm, previous work by morozov and siggia  <cit>  demonstrated that cumulative contact pwms could be used as informative priors in the task of scanning genomic sequences.
in this work we found that our algorithm outperforms the aforementioned method for homology models of tfs and displays a comparable performance with crystallographic structures.
this fact gives relevance to the statistical models that can be generated with dnaprot and preliminary tests while scanning genome-sized sequence sets with the model built for purr confirm this expectation.
overall, our study adds new insights to the challenging problem of estimating protein-dna binding specificities from structural complexes alone.

 the specific recognition of genomic cis-regulatory elements by nucleic acid binding proteins is of critical importance for many vital processes such as dna replication and repair, mrna translation and transcriptional regulation. specific tf-dna interactions, as well as protein-protein contacts established between closely located tfs at promoter regions, are key determinants of the concerted expression of genes in response to external stimuli. the task of determining the relationships among tfs and their targets has been widely addressed experimentally and computationally. of great relevance and use in recent years are computational approaches, based on the construction of statistical models that characterize the dna-binding preferences of tfs, which are in turn used to scan genomic sequences in an effort to identify new putative binding sites. the probabilistic models more commonly used in such computational approaches are position weight matrices obtained from multiple alignments of known binding sites. this approach is limited to tfs with a sufficient number of experimentally identified binding sites, for which reliable statistical models may be built. an alternative approach would be to predict dna operator sites which are compatible with the mode of binding of a given tf. this approximation would depend on two key components: i) the knowledge of the protein and dna residue positions involved in binding at the spatial level and ii) a method to evaluate the compatibility of different dna bases and amino acids to interact. the computational analysis of available protein-dna complexes has resulted in several papers describing the characteristics of the amino acid-base interactions that determine binding specificity, the different types of readout mechanisms involved in dna recognition and the evolutionary conservation of residues located at contact interfaces. direct readout is associated to recognition through contacts established between atoms from amino acid side-chains and nitrogen bases. indirect readout, on the other hand, is mediated by the contribution of residues of the proteins and dna which are not in direct contact and conformational changes undergone by dna upon protein binding. these reports show that both direct and indirect readouts significantly contribute to specific protein-dna recognition. with respect to direct readout, it is estimated that about two-thirds of all protein-dna interactions are van der waals contacts which do not generally confer sequence specificity, with the exception of the hydrophobic interactions involving the c atom of thymine. on the contrary, hydrogen bonds are the major source of specific interaction contacts. two-thirds of the hydrogen bonds formed between amino acids and bases are bi-dentate or complex interactions providing a great specificity. non-classical c-h ···o hydrogen mediated links have also been found at protein-dna contact interfaces, albeit their energetic contributions to recognition are not fully understood yet. these interactions are inherently weak individually, which means that cumulative effects are indispensable to significantly account for specificity, as occurs with hydrophobic links. water mediated bridges, though common in interaction interfaces, are mostly used as gap-fillers and are also engaged in buffering unfavorable electrostatic repulsions between interacting atoms at the interface. electrostatic interactions, most of which are established between the protein main chain and the sugar/phosphate backbone, do not significantly contribute to specific recognition, although they play important roles in the transition from the unspecific to the specific complex. to date just a few reports have been published aiming at finding putative tf target sequences using structural information. there have been some attempts to apply physical potential functions – i.e., in the form of atomic force fields – to estimate the energy of interaction and the relative contribution of direct and indirect readout mechanisms. however, these physical models do not appear to significantly outperform simpler statistical methods. examples of these methods are approaches oriented to extracting structural information from family-wise comparisons, building a statistical model of the interface derived from crystallographic structures and binding sites of proteins that belong to the same family. other procedures infer statistical potentials of interaction at the residue-base level from datasets of known complexes, which are used to examine the compatibility between the protein and its putative sites. however, there is some evidence suggesting that the performance of these methods may be limited by the simplicity of their interaction potentials, calculated at the level of cα atoms for protein residues. moreover, the increasing number of high-resolution structures for protein-dna complexes opens a new door towards a more detailed study of their contact interfaces. indeed a recent report confirms that atomic details improve the ability of physical potential functions to predict sequence-specific protein-dna interactions. in this work we constructed position weight matrices that capture the binding specificity of transcription factors, based on information extracted from the protein data bank. three atomic preference matrices, for hydrogen bonds, water-mediated hydrogen bonds and hydrophobic interactions, were derived from a non-redundant training set of complexes annotated in the pdb. these matrices were used to make explicit atomistic representations of hydrogen and hydrophobic bonds as well as the contribution of water molecules at interfaces, which gives us the opportunity to score the direct readout. this contribution is combined with empirical estimations of dna deformation, in order to calculate a potential of binding that includes both direct and indirect readout contributions. we evaluated the performance of our algorithm in a set of bacterial and eukaryotic tfs which have been co-crystallized bound to dna, and in most cases our results proved to be as good as or better than those obtained with the structure-based cumulative contact method by morozov and siggia. in addition, two tf homology models were analyzed in detail and used to predict their dna binding motif, after sampling side-chain rotamers at their contact interfaces. in this case the results we obtained were significantly better than those returned by the reference method, which indicates that our algorithm could be suitably used to study tfs of unknown structure starting from structural models. we also discuss the strengths and limitations of our approach that might potentially be used for tfs with few or no experimentally characterized binding sites.  protein-dna interface atomic contacts and interaction preferences starting from a non-redundant training set of crystallographic protein-dna complexes culled from the pdb, we constructed atomic frequency matrices for hydrogen bonds, water-mediated hydrogen bonds and hydrophobic contacts. a close inspection at the information embodied in these frequency matrices reveals some interesting features of protein-dna interfaces. for instance, as claimed in previous reports, we found that arginine is the major source of hydrogen bonds, with a marked preference to interact with guanine. however, in the present study we divided this preference in pairs of interacting atoms, finding that groups nh and nh of arginine establish 88% of the hydrogen bonds with atoms n and o from guanine. histidine also showed a marked preference towards guanine, with atom ne interacting in out of hydrogen bonds with atoms o and n from the nitrogen base. overall, we found interface hydrogen bonds in our training set. hydrogen bonds reported by hbplus between pairs of atoms in the side-chains of amino acids and bases in the set of protein-dna complexes extracted from pdb were computed. rows correspond to the hydrogen bond donor/acceptor atoms from the amino acids and columns to the hydrogen bond donor/acceptor atoms from bases with their corresponding pdb atom codes. each cell indicates the number of contacts found for a given pair of atoms. in contrast with hydrogen bonds, the landscape of hydrophobic interactions is quite different. as previously reported, the main discriminatory group regarding specific recognition through hydrophobic interactions is the methyl group of thymine. accordingly, the c group accounted for 20% of all hydrophobic interactions found in our training dataset, being the main source of contacts for all amino acids when compared with the other nitrogen bases. with the exception of c we found no obvious interaction preferences between amino acid and bases. this means that the signal to noise ratio of the calculated preferences is lower than the one observed for hydrogen bonds, and are therefore less informative. thus, in our subsequent analyses, we decided to take into account only the contribution of the methyl group of thymine during the evaluation of the interface between a tf and a given dna sequence. we also generated frequency matrices for water-mediated hydrogen bonds, finding a total of atomic interactions. the overall observed interaction propensities are quite similar to those in the hydrogen bond matrix, which is probably related to the gap filling role of water that supports the existence of long-distance hydrogen bonds. in addition, the buffering role of water allows the formation of electrostatically unfavorable hydrogen bonds, adding new interaction propensities absent in the hydrogen bond matrix. for instance, we found eight water-mediated hydrogen bonds between arginine side-chain nitrogen groups and cytosine n a survey of the quality of the atomic interaction matrices the redundancy of the data in the pdb and the bias towards specific protein folds considerably tangle the efforts aiming at exploiting the wealth of this database, limiting the scope of the results. in order to overcome these problems, we decided to avoid redundancy as much as possible in our training set, taking care of withholding a sufficiently informative set of pdb entries. we also planned a thorough bootstrap assay to estimate the quality of the atomic preferences extracted from the training dataset in this bootstrap tryout we resampled with replacement a thousand subsets of entries, randomly excepting a considerable fraction of the entries included in the initial training set as described in methods. by leaving out almost a third of the total entries, we measured the degree of overtraining due to residual redundancy in the initial set, as well as the statistical relevance of the interaction preferences reckoned in the matrix building process. a straightforward way of addressing this issue is by finding out whether the matrices constructed using the truncated training sets were good enough to suitably evaluate the direct readout for entries used to build the corresponding bootstrap matrices, though not failing to score the crystallographic complexes that were excluded in the resampling step. as may be seen in figure the exclusion of 30% of the entries in the training set does not significantly affect the discrimination capability of the matrices. here we assess the scores obtained when evaluating the bootstrapped and excluded datasets with the bootstrap matrices – in the abscissa -, against the generic or shuffled matrices. we found in the scatter plots and the regression analysis that the scores obtained for the assay of the bootstrapped and excluded datasets with the bootstrap matrices correlate very well – i.e., r = and r = respectively -, with the values obtained when assaying those same datasets with the generic matrices. as anticipated, the correlation entirely disappears when the analysis is done with the shuffled matrices – two bottom charts of figure direct and indirect readout mechanisms in protein-dna complexes the information embodied in the atomic interaction matrices presented in the previous section can be used to estimate the contribution of direct readout in protein-dna recognition. however, it has been shown that a mechanism of indirect readout also plays a relevant role in this process, at least for some tfs, and was therefore considered in this work. in particular, we modeled indirect readout as the cost of threading a nucleotide sequence into a fixed dna backbone. in order to integrate both direct and indirect recognition mechanisms, we designed a saturating mutation strategy, which is the kernel of the dnaprot algorithm. briefly, the algorithm iteratively evaluates the interaction potential of a given tf as the docked nucleotide sequence mutates in a 4n space, for a dna duplex of length n. this renders a structure-based position weight matrix that can be used to scan genomic sequences. a full description is provided in the methods section, but it is important to note that direct and indirect readout scores are linearly combined by means of a deformation weight, d. as d gets larger, the relevance of indirect readout increases. a study of the performance of our method is shown in figure where it can be observed that the indirect readout contribution to the ability of correctly scoring cognate sites of crp and narl is critical, since the predictive competence increases with large d values, as revealed by the increasing areas under the curve in the roc plots. these outcomes agreed with previous reports in the literature claiming the central role of deformation energy in the site recognition mechanism of crp and narl, which bend dna 90° and 42° respectively. purr and dnaa display an opposite trend, since their largest aucs are obtained with d values less than. while purr is known to bend dna with an angle of 45°, the reported bending angle of dnaa is just 28°. this means that although in the first two cases the contribution of deformation can be considered essential, in the cases of purr and dnaa optimal predictive values are obtained with intermediate values of d, also including an important contribution of direct readout. the subtle changes observed for the auc in the roc plots of figure in fact correspond to significant variations in the classification performance of our method. this can be seen in the results reported in table for the scoring of purr known binding sites. the comparison of the pwms obtained with d values of and, those for which we obtained the best and worst performance in figure 3c, reveals that setting d to its optimum value causes the p-values of cognate sites to be reduced by one or two orders of magnitude. this improvement corresponds to a reduction of 103- false positives extracted along with true sites in a search set of 106bp. it is worth noting that, using the optimum d value in the construction of the model for this tf and setting the sensitivity to 80%, the number of false positive sites recovered with the true sites is as low as in a genome-sized search set. the known binding sites of purr reported in regulondb were scored using dnaprot, maintaining w fixed to and setting d to and. these values of d correspond to the parameterization that renders the worst and the best pwms for this transcription factor respectively. the first column corresponds to the binding sites, the second and fourth columns report the scoring obtained with dnaprot and the third and fifth columns show the p-values of the scores obtained in a site recovery assay in a background set of random sequences. scoring crystallographic protein-dna complexes with the dnaprot algorithm in this section we present the results of scoring a diverse collection of protein-dna complexes solved by x-ray crystallography. the test set includes eight transcription factors bound to cognate dna sequences, of which four are prokaryotic and the other four eukaryotic, belonging to eight different superfamilies according to scop. in this experiment we took the coordinates of each of these interfaces and compute structural position weight matrices that approximate their binding specificity. two types of pwms are calculated here: i) readout matrices, derived using the dnaprot algorithm outlined in figure and ii) cumulative contact matrices, obtained by simply counting the contacts between protein side-chains and nitrogen bases, as previously described by morozov and siggia. these pwms can then be compared with cognate matrices built from experimentally determined binding sites, by means of local alignments. further details are provided in the methods section, but it is relevant to note that readout pwms are computed using knowledge-based atomic potentials, whereas cumulative contact pwms are calculated with the assumption that the consensus dna sequence is part of the pdb complex. figure includes the sequence logos derived from the computed pwms; the expectation values of pwm comparisons are shown in table together with details of the test set. scop v in the first column we include the name of the structural superfamily of each tf according to scop. in the second column the name of each tf is followed by the pdb identifier of the structure used enclosed in brackets. the following two columns contain the resolution and the r-value of the crystals respectively. the last two columns report the e-values resulting from stamp local alignments between the structure-based pwms obtained with the method of morozov and siggia or our approach and the corresponding cognate pwms, extracted from the literature. parenthesized e-values were calculated after sampling interface rotamers with dnaprot. the pwms generated by dnaprot starting from crystallographic structures are comparable or better than those produced with the methodology of cumulative contacts by morozov and siggia. this is clear in the results of table where, with the exception of purr and leu the statistical significance of our readout matrices is higher or of the same order as that of the contact matrices. an inspection of figure also shows that for the three prokaryotic regulatory proteins, crp, narl and metj, the sequence logos obtained by both methods are quite similar. however, except for the case of metj, for which the significance of the pwms generated by both methods is of the same order of magnitude, for crp and narl the matrices generated by our algorithm are two orders of magnitude better. on the contrary, although the logo of purr is correct, the e-value it gets when compared with the expected logo is several orders of magnitude worse, as a consequence of failing to identify consensus nucleotides. the results for our selection of eukaryotic tfs show close similarities for the statistical significance of the pwms generated by both methods. however, for rap sampling of interface rotamers caused an improvement of the readout matrix that outperforms the contact matrix by one order of magnitude,. something similar happened for the matrix generated for purr, though in this case the improvement of three orders of magnitude of the significance is still worse than the value of the contacts pwm. footprinting comparative models of protein-dna complexes with dnaprot in this section we further test the ability of dnaprot readout matrices to capture the binding specificity of tf-structural models. here we generated homology models for two transcription factors: fnr, a global regulatory protein in e. coli, and giant, a regulatory protein involved in the early development of d. melanogaster. as in the previous section, we also computed cumulative contact pwms for comparison purposes, and we evaluated both structure-based matrices by aligning them to cognate matrices built from operator sites reported in the literature. we now present the results of both models in more detail. the protein sequence of fnr can be confidently aligned to the sequence of e. coli crp with tfmodeller, despite an overall sequence identity of just 22%, providing a reliable frame for comparative modeling. however, several residue substitutions occur at the binding interface, suggesting that the readout might be changing from the template to the model. indeed we found that four out of ten interface residues are mutated in the dimeric model, as annotated in table for this reason we sampled different rotamers for the side-chains of the interface residues and selected those that yielded best readout scores using the dnaprot algorithm. as shown in figure this limited sampling was enough to improve the readout pwm, increasing the number of correct consensus nucleotides in the sequence logo from six to nine. this improvement of the model is also related to an increase of the statistical significance of the readout matrix of one order of magnitude with respect to the cumulative contact matrix, as shown in table a further improvement is also possible by deriving a contact matrix from the refined homology model, which causes an additional increase of the e-value to e- the multiple alignment, figure 5d, also suggests that the interface differences inferred from the comparative model correlate with the sequence conservation of the dna-binding domain among fnr and crp homologous sequences. scop v winged helix leucine zipper domain in the first column we include the name of the structural superfamily of each tf according to scop. the name of each model is followed by the pdb identifier of the template used to guide model building, enclosed in brackets. the following three columns contain the ratio of contact interface identity, the resolution and the r-value of the crystals respectively. the last two columns report the e-values resulting from stamp local alignments between the structure-based pwms obtained with the method of morozov and siggia or our approach and the corresponding cognate pwms, extracted from the literature. parenthesized e-values were calculated after sampling interface rotamers while bracketed e-values correspond to cumulative contact pwms computed over the refined dnaprot complex. the second study of a homology model was performed with giant, an eukaryotic leucine zipper transcription factor that has been studied previously. this protein can be modeled based on human cebpb, with a sequence alignment that covers 92% of the dna-binding domain. although the overall sequence identity is only 13%, six out of interface residues, as labeled by tfmodeller, are conserved. the first evaluated giant model yielded a pwm that captures only four correct consensus positions, as seen in figure as in the case of fnr, we refined the binding interface of giant by sampling rotamers and selecting those with best readout scores. the refined model yields a structure-based pwm that includes eight consensus nucleotides, improving the preliminary model results in two orders of magnitude in terms of e-value. this sampled model also suggests a different set of interface residues responsible for dna sequence discrimination. the figure shows an alignment of giant and cebpb homologous sequences that seem to cluster in two different modes of binding to dna. relative contribution of atomic interactions at dna recognition interfaces the dnaprot algorithm models the direct readout mechanism by labeling hydrogen bonds, hydrophobic interactions and water-mediated hydrogen bonds at the interface, which are evaluated by means of log-likelihood preference matrices. therefore, the data generated for the eight tfs presented in table can be further analyzed with the aim of calculating the relative contribution of each of these atomic interactions to dna recognition. as shown in table we found four cases in which recognition is mediated exclusively by hydrogen bonds, while the remaining tfs seem to employ two or more types of atomic interactions to drive dna binding. of these, perhaps the most interesting cases are narl, which has a large hydrophobic contribution, and zif with 30% of the interface score contributed by water-mediated hydrogen bonds. repeating these calculations on the set of non-redundant complexes, we found that 78% of dnaprot interface scores are contributed by hydrogen bonds, while water-mediated hydrogen bonds account for 16% and hydrophobic interactions are responsible for 6% of the total. table 5b shows a more detailed analysis of narl and zif in which structure-based pwms derived from only one type of atomic interactions are compared with the same cognate pwms in table these last results show that the combination of explicit representations of the relative contribution of the different types of interaction improves the predictive competence of dnaprot. a) transcription factors in table where further analyzed by assessing the fraction of the interface binding score, obtained after summing all log-likelihood atomic scores, contributed by each interaction type. the number of atomic interactions found by the dnaprot protocol in the native pdb structure is indicated in the rightmost column. b) a further analysis of two notable examples, narl with an important hydrophobic contribution and zif with a remarkable contribution of water-mediated hydrogen bonds, in which structurally derived pwms that consider exclusively h-bonds, hydrophobic contacts or water-mediated h-bonds are compared with default dnaprot pwms, which integrate all three atomic interactions. discussion the field of computational prediction of regulatory sites in genomic sequences has been dominated by methodologies relying on prior knowledge of experimentally characterized binding sites. only a limited group of studies has intended to generate structural representations of the sequence recognition interface, using approaches that somehow approximate the physical processes involved in dna recognition. some of these studies have dissected binding interfaces at the residue level, while others use atomic-based descriptions. though some efforts employ force fields for the calculations, most of these methodologies rely on pre-compiled interaction preferences. moreover, there are also conceptually simpler approaches that rely on counting contacts at the interface that generally work under the assumption that the nucleotide sequence captured in a crystallographic protein-dna complex is the consensus. these contact-based approaches, particularly the cumulative-contact method of morozov and siggia, promise to be successful in binding site searches at the genomic scale if the structure of a bound tf is at hand. extraction of protein-dna atomic interaction preferences from the pdb the dnaprot algorithm introduced here incorporates an atom-based representation of protein-dna interfaces and explicitly integrates both direct and indirect readout. though the thermodynamic contribution of interface desolvation and conformational entropy of interacting groups have been considered in some approaches using complex formulations, our algorithm, outlined in figure is concerned with binding specificity and does not account for those contributions, like other related reports. however, our methodology attempts to approximate these phenomena by including explicit representations of binding interface water molecules and side-chain rotamer sampling of interacting amino acids, which have important implications in our results. the hydrogen bond and hydrophobic interaction matrices used by dnaprot were derived from a non-redundant set of complexes, naturally surpassing the training sets used in preceding studies, as more structural data is now available at the pdb. an important assumption behind the construction of this data set is that rules governing specific recognition of dna by proteins are generally due to the conformational restrictions imposed by the double helix. therefore, we chose to collect a comprehensive training set instead of using a superfamily-specific or a tf-exclusive set, in which the scantiness of information would have weakened the derived statistical models. in addition, we built an atomic contact matrix that explicitly accounts for water-mediated hydrogen bonds in protein-dna interfaces, which constitutes a novel contribution of this work. the data obtained in the construction of our contact propensity matrices, such as the hydrogen bond matrix in table constitute an update of the pioneer work by mandel-gutfreund et al. and luscombe et al. as we found similar trends in the atomic interaction propensities in our enlarged dataset. in addition, this work presents quantitative arguments, which demonstrate the importance of including hydrophobic and water-mediated hydrogen bonds in the list of interactions that contribute to specific recognition of nucleotide sequences. assessment of the quality of the atomic interaction matrices keeping in mind the problem of data bias in the pdb, we planned an in-depth bootstrap tryout to assess the reliability of the information contained in our contact matrices. the exclusion of almost a third of the entries in the original training set yielded atomic interaction matrices with only a slightly reduced ability to correctly evaluate infrequent atomic contacts. these results confirm that the preference matrices derived in this work are biologically significant and that the atomic associations found are not significantly affected by overtraining nor by redundancy, since bootstrap matrices scored the excluded training sets almost as well as the bootstrapped ones. on the contrary, the results for the negative control obtained while repeating this analysis with shuffled matrices, show completely scattered distributions with no correlation. these outcomes support the consistency of the interaction preferences extracted from the training set used to build our matrices. they also corroborate that the inclusion of new protein-dna complexes, regularly deposited in the pdb, will not cause considerable variations of the interaction preferences, except for very unusual atomic contacts. evaluation of direct/indirect readout contributions to protein-dna recognition the rationale for using atomic potentials is that many interface contacts are either complex or bidentate, and therefore cannot be properly accounted for by residue-based approaches. in those cases the contact pairing ratio is not 1: as the same amino acid might be in contact with two base steps simultaneously or binding a given nitrogen base through multiple groups. these drawbacks could be resolved by using explicit atomistic representations of the mode of binding, which have been the goal of a group of papers in this field. however, using atomic potentials requires high-quality atomic structures of protein complexes. in addition, atomic-detail approaches such as dnaprot might be more affected by cases in which side-chains rearrange upon mutation of the bound nucleotides, which is the reason that it might be necessary to sample interface rotamers while doing comparative modeling, as explained in the results section. further sampling, including backbone movements, might be required in certain cases, as already envisaged by havranek et al.. in fact, previous works have already recognized that homologous protein-dna interfaces change their docking geometry as their sequences diverge. besides the contribution of the direct readout to sequence recognition, we also considered in our model a mechanism of indirect readout, as sequence-specific dna deformation has been identified as key to dna recognition for many transcription factors. the algorithm presented in this work follows previous efforts that model indirect readout as the cost of deforming a dna duplex. both readout mechanisms, shown to be critical for specific recognition in figure are linearly integrated into a single binding potential. the weights of both direct and indirect readout can be tuned for different transcription factors according to their docking mode, as the examples depicted in figure imply. this observation suggests that the performance of the dnaprot protocol can be improved using previous biological knowledge and this is certainly a desirable property. this data also insinuates that each tf might have its own balance between direct and indirect readout, although we cannot exclude the possibility that this value is a property of its structural superfamily. moreover, data in table gives further insight into the impact that correct weighting could have on genomic scale tf-binding-site prediction assays, as the number of false positives might be considerably reduced in order to obtain a more reliable set of predictions. assaying the predictive potential of dnaprot in crystallographic structures and homology models the predictive power of dnaprot, evaluated in tables and and figures and suggests that our readout pwms are, with one exception, comparable to or better than those generated by a well established reference methodology, the cumulative contact pwm proposed by morozov and siggia and related to previous reports. it is worth remembering the reader that, unlike the reference method, dnaprot does not assume that the nucleotide sequence in the input pdb complex is the consensus; rather, it performs an in silico mutagenesis assay and evaluates 4n sequences using a readout function. as already mentioned, the conceptually simpler cumulative contact approach does not consider the contribution of indirect readout. this might explain the better results obtained with our method for crp and narl shown in table two tfs known to have an important contribution of indirect readout. in contrast, the example of purr shows that the dnaprot formulation of hydrogen bonds and hydrophobic interactions does not fully capture the array of atomic contacts at the interface, as the obtained readout sequence logo misses consensus positions that might correspond to a different type of interaction. despite this fact, the predicted pwm correctly identifies the remaining consensus positions except one, and the statistical significances obtained in the site recovery assay of the cognate binding sites of this tf depicted in table prove the predictive competence of the matrix built. notwithstanding the promising results obtained for crystallographic structures, the most valuable application of the methodology presented in this paper is found in the exercise of comparative modeling. two examples were modeled here, a well-characterized prokaryotic regulator and an eukaryotic factor previously studied in a related article. both examples demonstrate that rotamer sampling at the interface is necessary for obtaining optimal results, and we found that the readout formulation presented here can be effective in selecting the best rotamers. as shown in table and figure our methodology yielded better results than those obtained with the cumulative contact method for two homology models. the reason for this could be that comparative models usually contain errors in the assigned position of protein side-chains and, most importantly, they do not necessarily contain the consensus dna sequence. a combination of in silico saturating dna mutagenesis and interface side-chain sampling allows the dnaprot algorithm to partially overcome these problems. not only is sampling positive in modeling tasks, but we also found that the crystallographic structures of purr and rap yielded better readout pwms after resampling their interface side-chain rotamers. although the quality of the structures used by our methodology is of primary relevance, we could not find a clear correlation between r-factor, resolution or experimental techniques for entries in the pdb and the outcomes of our procedure. nevertheless, we often encountered protein-dna complexes where only a few hydrogen bonds or hydrophobic interactions can be identified using standard bond geometries, yielding only partial dna motifs. in such cases we found that the cumulative contact method seems to be less sensitive to the quality of structural data. further work needs to be done to explore whether side-chain sampling, and even backbone sampling, can help circumvent the limitations that data quality imposes on the performance of readout pwms. evaluation of the contribution of different interaction types to the recognition process the atomistic foundation of our methodology also gives us the possibility of exploring the relative contribution of each interaction type to the dna recognition process. the analysis of a non-redundant set of more than complexes proves that water-mediated hydrogen bonds are the second source of specific interactions at contact interfaces, which raises a warning regarding the usual exclusion of these interactions in structural studies trying to model protein-dna interaction interfaces. the exclusion of water overlooks a wide group of highly informative contacts, mainly long-distance and unfavorable hydrogen bonds that constitute novel pairings absent from the hydrogen bond preference matrix. the example of zif analyzed in table 5b, shows that an important part of the information content of the pwm obtained for this tf corresponds to water-mediated interactions, and the inclusion of this information considerably increases the statistical significance of the readout matrix. something similar occurs for hydrophobic interactions; despite being relatively infrequent sometimes their contribution to binding might be central. this is the case of the matrix obtained for narl, in which hydrophobic links account for most of the information content of the structurally-derived pwm, despite being the least frequent interface atomic interaction. in this last example, without a reliable inclusion of those infrequent interactions, our model for this tf would have been worthless.  in summary, our results suggest that the dnaprot algorithm, together with the set of atomic interaction matrices obtained in this work, have useful applications. the matrices contain biologically meaningful information that confirm and enrich previous reports at the atomic level of interaction. in addition to the uses presented in this paper, these matrices could be taken to estimate specificity of binding or as a guide for crystallographic studies of protein-dna complexes. with respect to the algorithm, previous work by morozov and siggia demonstrated that cumulative contact pwms could be used as informative priors in the task of scanning genomic sequences. in this work we found that our algorithm outperforms the aforementioned method for homology models of tfs and displays a comparable performance with crystallographic structures. this fact gives relevance to the statistical models that can be generated with dnaprot and preliminary tests while scanning genome-sized sequence sets with the model built for purr confirm this expectation. overall, our study adds new insights to the challenging problem of estimating protein-dna binding specificities from structural complexes alone. 
we first measure the ability of our method to correctly estimate the binding specificities of eight prokaryotic and eukaryotic tfs that belong to different structural superfamilies. secondly, the method is applied to two homology models, finding that sampling of interface side-chain rotamers remarkably improves the results. the specific recognition of genomic cis-regulatory elements by transcription factors plays an essential role in the regulation of coordinated gene expression. our algorithm combines information regarding the direct and indirect readout of dna into an atomistic statistical model, which is used to estimate the interaction potential. studying the mechanisms determining binding specificity in protein-dna interactions is thus an important goal. our results demonstrate that atomic-detail structural information can be feasibly used to predict tf binding sites. thirdly, the algorithm is compared with a reference structural method based on contact counts, obtaining comparable predictions for the experimental complexes and more accurate sequence motifs for the homology models. most current approaches for modeling tf specific recognition rely on the knowledge of large sets of cognate target sites and consider only the information contained in their primary sequence. the computational method presented here is universal and might be applied to other systems involving protein-dna recognition.
 the specific recognition of genomic cis-regulatory elements by transcription factors plays an essential role in the regulation of coordinated gene expression. studying the mechanisms determining binding specificity in protein-dna interactions is thus an important goal. most current approaches for modeling tf specific recognition rely on the knowledge of large sets of cognate target sites and consider only the information contained in their primary sequence.  here we describe a structure-based methodology for predicting sequence motifs starting from the coordinates of a tf-dna complex. our algorithm combines information regarding the direct and indirect readout of dna into an atomistic statistical model, which is used to estimate the interaction potential. we first measure the ability of our method to correctly estimate the binding specificities of eight prokaryotic and eukaryotic tfs that belong to different structural superfamilies. secondly, the method is applied to two homology models, finding that sampling of interface side-chain rotamers remarkably improves the results. thirdly, the algorithm is compared with a reference structural method based on contact counts, obtaining comparable predictions for the experimental complexes and more accurate sequence motifs for the homology models.  our results demonstrate that atomic-detail structural information can be feasibly used to predict tf binding sites. the computational method presented here is universal and might be applied to other systems involving protein-dna recognition. 