Abstract
Background
A computational method (called p53HMM) is presented that utilizes Profile Hidden Markov Models (PHMMs) to estimate the relative binding affinities of putative p53 response elements (REs), both p53 singlesites and clustersites. These models incorporate a novel "Corresponded BaumWelch" training algorithm that provides increased predictive power by exploiting the redundancy of information found in the repeated, palindromic p53binding motif. The predictive accuracy of these new models are compared against other predictive models, including position specific score matrices (PSSMs, or weight matrices). We also present a new dynamic acceptance threshold, dependent upon a putative binding site's distance from the Transcription Start Site (TSS) and its estimated binding affinity. This new criteria for classifying putative p53binding sites increases predictive accuracy by reducing the false positive rate.
Results
Training a Profile Hidden Markov Model with corresponding positions matching a combinedpalindromic p53binding motif creates the best p53RE predictive model. The p53HMM algorithm is available online: http://tools.csb.ias.edu webcite
Conclusion
Using Profile Hidden Markov Models with training methods that exploit the redundant information of the homotetramer p53 binding site provides better predictive models than weight matrices (PSSMs). These methods may also boost performance when applied to other transcription factor binding sites.
Background
The p53 protein plays a crucial role in cancer suppression in the human body. In response to cancerinducing, DNAdamaging stress conditions, the tetrameric p53 proteins can activate different pathways that lead to DNA repair, cell cycle arrest, inhibition of angiogenesis, and apoptosis [1]. A highly degenerative, palindromic consensus DNA binding site, consisting of a halfsite RRRCWWGYYY, followed by a variable length spacer, then followed (almost always) by a second halfsite RRRCWWGYYY sequence, has been discovered for the protein, where R is a purine, Y a pyrimidine, W is either A or T (adenine or thymine) and G is guanine and C is cytosine (see Figure 1) [2,3]. By labeling each quartersite RRRCW as → and its reversecomplement WGYYY as ←, the first discovered p53 consensus sequence can be graphically represented by → ← spacer → ←. This configuration of the four quartersites is often referred to as the headtohead (HH) orientation, and represents the vast majority of experimentallyvalidated p53 binding sites to date.
Figure 1. Original Data from ElDeiry et al., Used To Define The p53 Consensus Binding Site. The original DNA fragments collected from a genomewide, p53antibody immunoprecipitation, that were used to define the headtohead (HH) p53 Consensus Binding Site, are graphically presented [3]. The yellow columns corresponding to the 1^{st }and 2^{nd }halfsites were used to define the consensus p53 motif. The p53 binding site is highly degenerative. Within the yellow columns, notice that 7 of the 20 DNA target sites (35%) had at least one nucleotide insertion (green), deletion (red), or both (magenta) relative to the discovered 10 bpspacer10 bp consensus. Since insertions and deletions throw off the reading frame of a weight matrix, any PSSM approach will inherently misscore at least 35% of these 20 sites. Alignments of the 160 experimentally validated p53 binding sites also reveal that any PSSM approach would inherently misscore at least 30% of them as well. Another observation is that additional p53 halfsites are immediately adjacent (in yellow) to the ones used to define the consensus in 15 of the 20 target sites (75%). Since the genomewide immunoprecipitation study was designed to pull down the highest affinity sites, the fact that 75% of the target sites are actually p53 clustersites is the first indication that clustersites of 3 or more halfsites confer higher binding affinity [22].
The degeneracy of the p53RE
In the influential paper "Definition of a Consensus binding Site for p53", by ElDeiry et al., 7 of the 20 DNA target sites (35%) used to form the headtohead (HH) p53 consensus sequence had at least one nucleotide insertion or deletion relative to the discovered 20 bp consensus after proper alignment (see Figure 1) [3]. Alignments of the roughly 160 experimentallyvalidated p53 binding sites to date also show that approximately 30% of presently known sites have at least one nucleotide insertion or deletion relative to the consensus matrix [4]. Discovery of p53 binding sites with such degeneracy cannot be reliably made with a PSSM approach, since prevalent insertions and deletions in the consensus sequence misalign the PSSM reading frame, and lead to improper scoring. Therefore, PSSM binding site discovery algorithms inherently misscore at least 30% of the known p53 binding sites.
PHMMs can model nucleotide insertions and deletions
Profile Hidden Markov Models provide a coherent theory for probabilistic modeling of degenerate binding sites where random nucleotide insertions into and deletions from the motif are tolerated at certain positions [5,6]. Natural selection suggests that critical nucleotides are conserved over evolutionary time, while noncritical nucleotides (including tolerated insertions in the motif) are not conserved. The matchstate emissions of the PHMM serve to model the critical positions in the motif with their observed nucleotide frequencies. The additional hidden deletion and insertion states at each position enable the model to train for (relatively rare) observed deletions and insertions at different positions in the motif (see Figure 2). Although the probability of any particular insertion or deletion of a nucleotide at a certain position in a functional motif may be rare, the accumulated probability over all the positions in the motif that an insertion or deletion event may occur can be significant. The training set of observed insertions and deletions serves to finetune the model to be properly sensitive to tolerated deviations from the most prevalent consensus motif. The main strength of the PHMM is this trained flexibility to properly model variable length motifs. The major drawback is that more data is required to train the extra parameters not found in weight matrices (PSSMs).
Figure 2. The Topologies of p53 Singlesite and Clustersite Models. (a) A Profile Hidden Markov Model (PHMM) contains three hidden states for each position in a sequence motif of length n: a match state (green squares), an insertion state (orange diamonds), and a delete state (gray circles). The arrows represent allowed transitions between states and have associated probabilities. The match and insertion states also have associated nucleotide emission probabilities. The first and last insertion states (I0 and In) and associated transitions (in red) are shown for completeness. However, they are not present in the p53 models since they are replaced by FIM and FEM models. (b) The topology of the Finite Emission Module (FEM) of length N allows the ability to model any distribution of spacerlengths between 1 and N. For the p53 models, the model and background probabilities within the FEM modules are identically uniform so that there is nocost for spacerlengths between 1 and N, and are referred to as "nocost FEMs". (c) The topology of the Free Insertion Module (FIM) allows for the ability to model an exponentially decaying distribution of spacerlengths. However, by setting the model and background probabilities to identically uniform, the FIM can model any sequence of infinite length with no associated cost to the overall score (hence the word "Free"). (d) The main components of the p53 singlesite model are the left and right halfsite PHMMs, which potentially contain corresponding positions between them. These two halfsite models are separated by a nocost FEM model that limits the length of any intervening spacer sequence to 20 bp. The halfsite models are also wrapped by two FIMs that allow the Viterbi algorithm to find the best matching motifs anywhere in the candidate sequences. (e) The topology of the p53 clustersite model consists of a single PHMM that models a general halfsite, and two backtransitions that allow for modeling an infinite number of halfsites within the clustersite. The backtransition through the nocost FEM14 model limits the spacersequence between the halfsites to lengths ≤ 14 bp.
Using PHMMs to estimate binding affinities
Like weight matrices (PSSMs), Profile Hidden Markov Models can be used to estimate the relative binding affinity of a protein for a particular binding site sequence [7]. Under ideal conditions, the logodds scores G^{s}(x) that a Profile Hidden Markov Model (trained on training set S) calculates for any candidate site x is directly proportional to the free energy ΔG(x) of the TFprotein binding to that candidate site [see Additional file 1 for details] [79]. The logodds scores are given by:
where we define:
With these definitions, and assuming independence of positions, we have:
The dynamic programming forward and backward algorithms are used to calculate the probabilities P_{hmm}(x) and P_{hmm}(j, b). These two probabilities are calculated by summing up the probability of observing the sequence x, and the base b at position j, for all the paths through the linear PHMM, respectively. The dynamic programming Viterbi algorithm is used to find the best alignment of the candidate site x to the bindingsite motif modeled by the PHMM. The best (optimal) alignment of the sequence x is obtained by finding the path through the PHMM that gives the highest logodds score for the sequence [8]. In the case of transcription factor binding sites, the logodds score of this optimal path (also called the Viterbi score) is commonly used to provide adequate approximations to the probabilities P_{hmm}(x) and P_{hmm}(j, b) [see Additional file 1 for details]. When using the Viterbi score for the probability P_{hmm}(x) we are assuming that there is generally only one major set of binding interactions between specific nucleotides and amino acids for a given proteinDNA complex, and that all other possible binding locations in the response element can be ignored.
Additional File 1. Supplementary Material. The Supplementary Material contains more theory behind modeling TFbinding sites with PHMMs, details of the Corresponded BaumWelch algorithm, and a proof that the PHMM logodds score of a TFbinding site estimates its relative binding affinity given certain assumptions.
Format: PDF Size: 201KB Download file
This file can be viewed with: Adobe Acrobat Reader
Training a PHMM with validated binding sites
Before a PHMM can be used to estimate the relative binding affinity for any putative binding site, the PHMM must be trained to properly model a functional binding site of interest. When training a PHMM for a particular motif, the goal is to choose the parameters of the model in order to maximize the likelihood of the sequences in the training set, without overfitting. Again, under ideal conditions the logodds score (loglikelihood ratio) G^{s}(x) to be maximized for the collection of binding sites in the training set is proportional to the estimated binding free energy ΔG(x) of these binding sites. When the state paths for the training sequences are not known, no known closed form solution exists for the parameter estimations [8]. The BaumWelch algorithm is the most commonly used iterative Expectation Maximization (EM) method to train the parameters of the model. The BaumWelch algorithm always climbs the gradient (to increase the combined scores of the training set) and uses the optimized dynamic programming forward and backward algorithms [8].
Results and discussion
A novel training method that boosts predictive power
To increase the predictive power of our p53motif PHMMs, we attempt to exploit the a priori knowledge that when proteins bind as homodimers or homotetramers, their corresponding binding sites typically have a palindromic, repeat, and/or reverse complement structure (see Figure 3). This prior knowledge can be used to correspond (fully or partially tie) the parameters between positions in order to exploit the inherent redundancy in the information of the motif. Within a set of corresponding positions, the updating of emission and transition probabilities can borrow strength from each other by sharing information. In addition, the degree of sharing of information for any set of corresponding positions can be optimized during training. The process of corresponding parameters can greatly reduce the parameter searchspace during the training of the model, and provide the ability to train for rare occurrence insertion and deletion events. This general technique has been effectively used when HMMs have been applied to speech and handwriting recognition problems, and has been referred to as parameter tying [10]. We introduce an extension to this method that allows for the setting or training for an optimal level of partial or full parameter tying. In the domain of proteinDNA binding sites, even if a palindromic, repeat, or reverse complement structure of a binding site is not known a priori, all the known structural motifs can be tested, and the structure can be discovered (inferred) from the ROC curve that maximizes predictive accuracy. For example, of the six structural models tested for the p53 binding motif, the combinedpalindromic motif that completely corresponds the four quartersites is the discovered motif, since it is the best classifier (see Figure 4).
Figure 3. The Four p53 Correspondence Motifs. The four correspondence motifs for the repeated, palindromic p53 RE are graphically represented. In the top three motifs, each line corresponds 2 synonymous positions. In the bottom motif, the previously independent halfsites are made corresponding (tied) by the yellow connecting lines so that now 4 synonymous positions are corresponded. The completely untied motif (not shown) has no correspondence, and thus no connecting lines, between any of the positions in the motif. (R = A or G, W = A or T, and Y = C or T. Position ã has the complement nucleotide emission distribution of a.)
Figure 4. Cross Validation with Receiver Operating Characteristic (ROC) curves reveals increased predictive power over weight matrices. 1000 iterations of 10fold randomsplit cross validation reveal that the most predictive models utilize the correspondence structures. The combinedpalindromic model is the best model since it contains roughly half as many parameters as the other three correspondence models. The positive set contains 160 experimentally validated p53 binding sites, and the negative set contains 40 bp random samples from the mononucleotide content of the training set. The true positive and false positive rates are calculated and plotted for all possible threshold values for each model. The predictive measure for comparing the curves is the AUC (Area Under the Curve). In all the PHMM models the insertstate emissions are fixed to the A, G, C, T nucleotide distribution of the training set. The best classifier uses the combinedpalindromic training motif. (Position ã has the complement nucleotide emission distribution of a).
The Corresponded BaumWelch algorithm
In order to include the prior knowledge of the structural motif (or in an attempt to discover it), a novel "Corresponded BaumWelch" algorithm is proposed to enforce or learn the optimal correspondence between expectations of parameters for corresponding positions after each iteration of the BaumWelch algorithm (see Methods). For example, assume that we have prior knowledge that a transcription factor protein binds to the DNA in homodimer form, where each monomer interacts with 5 DNA base pairs. Then a corresponding palindromic motif for the nucleotide positions would be: 1 2 3 4 5 5 4 3 2 1, while a reversecomplement palindromic motif would be: 1 2 3 4 5 (where ã has the complement nucleotide emission distribution of a). All the emission distributions for each of the five sets of synonymous positions would be made corresponding, as well as all the transition probabilities between synonymous positions. In this example, if all the parameters between synonymous positions were fully corresponding (tied), then the parameter search space would be roughly cut in half. The level of correspondence between the parameters for synonymous positions can be given a priori, or trained for if the training set is sufficiently large. One optimal level of correspondence, c, can be calculated for the whole motif (for all the corresponding positions), or a separate one can be found for each set of corresponding positions. (See Methods for details.)
Comparing the different p53 corresponding (structural) motifs
Since the 20 bptetrameric p53 binding site has a repeated and nested palindromic structure, different correspondence motifs can be constructed to train the PHMM models, and cross validation can be used to compare their predictive properties. The motifs that are compared are: the repeat or Tcoupled motif (1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10), the (reversecomplement) palindromic or Hcoupled motif (1 2 3 4 5 6 7 8 9 10 ), the independently (reversecomplement) palindromic or Qcoupled motif (1 2 3 4 5 6 7 8 9 10 , the repeated, fullypalindromic or combinedpalindromic motif (1 2 3 4 5 1 2 3 4 5 ), and the completely untied motif with no correspondence between any positions (see Figure 3) [11]. We perform 1000 iterations of tenfold randomsplit cross validation on each model to gain statistics on their predictive accuracy. The positive set contains 160 experimentally validated p53 binding sites from [4], and the negative set contains 40 bp random samples from the mononucleotide content of the training set. Then we utilize Receiver Operating Characteristic (ROC) curves in order to compare the predictive power of the classifiers in an unbiased, thresholdindependent (nonparametric) manner. This is achieved by calculating the true positive and false positive rates for all possible threshold values for each model. The summary statistic for comparing the ROC curves is the AUC (Area Under the Curve). AUC values lie somewhere between 1.0 and 0.5 (where an AUC of 1.0 would correspond to a perfect classifier, and an AUC of 0.5 would correspond to a classifier that is no better than random coin flipping.)
Training InsertState Emissions
A major consideration when training Profile Hidden Markov Models (PHMMs) is which parameters to train for at each position, and which parameters to fix at each position to the overall average. The more nonfixed parameters that must be trained for at each position in the motif, the more data that is needed to properly train the model. Ideally, a sufficiently large training set is available to be able to train for all the parameters in the PHMM at each position. Unfortunately, in the case of transcription factor binding sites, this is rarely the case. Typically, when using PHMMs to model DNA binding sites, both the insert probabilities and insert state nucleotide emissions probabilities are set to the binding site averages, since there are rarely enough examples of these rare occurrence events at a particular position to train those parameters for that position alone [12]. By corresponding (fully or partially tying) positions and in effect increasing the training data for each position, it may be possible to train the insertionstate emissions distributions for these corresponding positions. This could possibly boost predictive power of the models, if the p53 protein is selective as to which nucleotides can be inserted into the motif at certain positions without compromising the binding affinity of the site. A common example of such selective sequence insertions can be found in functional protein families, whereby hydrophobic or hydrophilic amino acid insertions may be tolerated at certain positions, provided that the insertions are present either in the core or at the surface of the protein, respectively, after folding. Notice that fixing the insertionstate emission distributions at every position to the aminoacid average for the whole sequence would be very inappropriate in this example.
The final results
The combinedpalindromic motif (1 2 3 4 5 1 2 3 4 5 ) performs on par with or better than all other structural motifs, although it contains comparably half the degrees of freedom (see Figures 4 and 5). In addition, all four of the structural motifs perform on par with each other. These results suggest that there exist correlations between the positions in the repeat, independently palindromic, and palindromic motifs, and that the combinedpalindromic motif leverages the correlations found in all of them. Furthermore, it can be seen that training the insertstate emissions per corresponding position also boosts the predictive power of all the models (see Figures 4 and 5). Analysis of the AUC measurements reveals some interesting features. Adding insertstate emission training to the base PHMM (with no motifcorresponded positions) has an AUC improvement of .923  .919 = .004, but with motif training has one of .937  .929 = .008. Adding motif training (motifcorresponded positions) to the PHMM when not insertstate emission training has an AUC improvement of .929  .919 = .010, but with insertstate emission training has one of .937  .923 = .014. Therefore the improvements are not additive. There is "positive synergy" when performing both motif training and insertstate emission training together that further boosts the predictive accuracy of the model. This observation confirms our hypothesis that training insertstate emissions can significantly boost the accuracy of the model after corresponding positions in the PHMM according to a bindingsite motif.
Figure 5. Cross Validation with Receiver Operating Characteristic (ROC) curves reveals increased predictive power when training insertstate emissions. All the PHMM models in this comparison train the insertstate emission distributions based on positional insertions occurring in the training set. Again, 1000 iterations of 10fold randomsplit cross validation reveal that the most predictive models utilize the correspondence structures. The positive set contains 160 experimentally validated p53 binding sites, and the negative set contains 40 bp random samples from the mononucleotide content of the training set. The true positive and false positive rates are calculated and plotted for all possible threshold values for each model. The predictive measure for comparing the curves is the AUC (Area Under the Curve). The AUC values improve for all the PHMM models compared to Figure 4, but not for the weightmatrix model (which does not use the insert states). The best classifier (with the combinedpalindromic training motif) was used for the p53HMM algorithm. (Position ã has the complement nucleotide emission distribution of a).
In addition, the more correspondence placed between the synonymous positions during each training iteration, the better the resulting classifier at that point in the training (results not shown). For this training set, all the combinedpalindromic models with fixed correspondence factors between c = 0.4 and c = 1.0 eventually converged to the same predictive model, although lower correspondence factors required more iterations to do so. All the models converged on correspondence factors between c = .98 and c = .999 when training for optimum correspondence. Therefore the best predictive model completely corresponds (ties) the four quartersites in a combinedpalindromic structure during each iteration of the training. Our published p53HMM algorithm is this best predictive model: trained on the dataset of 160 functional p53 REs, fully corresponding the data per position based on the combinedpalindromic structural motif, and training the insertstate emissions (see Figure 6).
Figure 6. The p53HMM Match and Insert Emissions. (a) The matchstate sequence logo for the combinedpalindromic p53 motif: 1 2 3 4 5 1 2 3 4 5 . (Motif position ã has the complement nucleotideemission distribution of a.) The height of each letter is made proportional to its frequency at each position, and the letters are sorted in descending frequency order. The height of the entire stack at each position is then adjusted to signify the information content (in bits) of that position [25]. The matchstate nucleotide positions 4, 7, 14, and 17 (motif positions 4, 7, , and respectively) are the most conserved and are the main points of contact with the p53 protein. (b) The insertstate sequence logo for the same combinedpalindromic p53model. These nucleotide insertions occur inbetween the nucleotide positions shown in part a. The specificity motif of the insertstate emissions is different from that of the matchstate emissions. (c) The HMM logo that combines parts a and b and state transition information into one graph. The wide, whitebackground stacks correspond to the match states in part a, while the narrow, redbackground stacks correspond to the insert states in part b. (A weakness of this HMM logo is that the insertstate stacks are so narrow that it is difficult to accurately see the stack specificity depicted in part b.) The yaxis is the same for all three graphs. However, the width of a stack in the HMM Logo is proportional to the expected contribution of that match or insert state to an emitted sequence of the model [26].
Validation of the p53HMM algorithm
The new p53HMM algorithm was used to screen for putative p53 binding sites in the endosomal compartment genes, which led to the discovery of a functional p53 site and a new p53regulated gene, CHMP4C [13]. The putative p53RE sequence AAACAAGCCC agtagcagcagctgctcc GAGCTTGCCC was predicted in the promoter region (497 to 460 bp) of the CHMP4C gene. The data from the chromatin immunoprecipitation and the luciferase reporter assays showed that p53 protein can bind to this sequence and induce CHMP4C gene expression. Additionally, analysis by p53HMM found an alternative putative p53 binding site in the LIF gene that corresponds to a 6 bp upstream shift of the downstream halfsite relative to the recently published putative site in intron 1 [14]. The p53HMM algorithm predicted the site GGACATGTCGGGACAGCTC, which matches the consensus RRRCWWGYYYRRRCWWGYYY perfectly except for the lowconserved position 10 and the gap ("", deletion) at position 16. A PSSM approach predicted the shifted site GGACATGTCGggacagCTCCCAGCTC, which is the best "gapless" p53 site in the region conferring p53 regulation, but it still matches the consensus very poorly with five mismatches (the putative spacer sequence is in lowercase) [14]. A few genes in the dataset of 160 functional p53 binding sites have a deletion relative to the consensus exactly between the wellconserved C and G as seen above, including the genes: EGFR, TYRP1, EEF1A1, HSP90AB1, and BAI1. This discovery of an alternative p53 binding site that better matches known functional sites, by modeling for observed insertions and deletions, highlights some of the advantages of the new p53HMM algorithm.
Special considerations for the p53HMM algorithm
Although the spacer within a p53 RE has been shown to greatly affect the binding affinity for p53 protein, the ability to properly quantify this effect for all possible spacers of lengths 0–21 base pairs has been elusive. Therefore like previous algorithms, we have chosen to initially ignore the spacers of the training set and putative REs [15]. We are able to ignore arbitrarylength spacers by inserting a nocost Free Insertion Module (FIM) between the two halfsites of the singlesite PHMM [16,17]. Similarly, we can ignore spacers with lengths between 1 and N base pairs by inserting a nocost Finite Emission Module (FEMN) between the two halfsites (see Figure 2). A prior p53 RE search algorithm (p53MH) was based upon a PSSM approach and a novel filtering matrix [15]. Unfortunately, the tables were not symmetric and the filtering table overfit the available data at the time. The combined result was that the p53MH method completely rejects 58 of the 160 experimentally validated sites to date (receiving a score of 0 out of 100, where 100 represented the maximum relative binding affinity). Additionally, some sites received very high scores approaching 100, while the reversecomplement received a score of 0, and viceversa. Due to these observations, we have purposely designed the p53HMM algorithm to be symmetric, so as to give identical scores for putative sites and their reverse complements. Secondly, we chose to abandon the filtering matrix to avoid overfitting the available data. A feature that we preserved from p53MH is the normalizing of scores by the highest possible affinity for the motif (×100), so that the highest possible normalized score is 100.
Modeling dependencies between positions
PSSMs assume that all nucleotide positions within the motif contribute independently to the binding affinity of the binding site, which has been shown experimentally to not always be the case [7]. Recent research has focused on modeling dependencies between positions in proteinDNA binding sites [18,19]. Typically Tree Bayesian Networks and Mixtures of trees have been used to attempt to model these dependencies between positions, which have been shown through cross validation to increase the predictive power of these models [18]. Our PHMM models do not attempt to model dependencies between the positions, however they can be extended to do so by using higherorder Profile Hidden Markov Models. Unfortunately, the ability to train for positional dependencies, and boost predictive power, is dependent upon the sampling size of the training set and requires larger training sets to train the extra parameters.
A novel p53 clustersite algorithm
Binding affinity measurements have been obtained for certain p53 clustersites of different lengths by mutating or truncating known p53 clustersites in the genes: DDB2, TP53i3, CKM, IGFBP3, and RGC (see Table 1 and Figure 7) [2023]. Based on the relative binding affinities of these p53 clustersites, we propose a new p53 clustersite algorithm that utilizes the trained PHMM to calculate and sum up the relative estimated bindingaffinities, above a certain threshold, of all viable fullsites in the cluster with a spacer of ≤ 14 bp or less (see Methods). This model predicts a linear increase in p53 binding affinity dependent upon the number of halfsites in the clustersite and the length of spacers between them. For example, for p53 clustersites with 2, 3, 4, 5, or 6 adjacent p53 halfsites, the number of possible fullsites with spacerlengths = 14 bp would be 1, 3, 5, 7, and 9, respectively. Let N be the number of halfsites in the clustersite, then the number of fullsites (to calculate binding affinities for and sum up) is given by the expression 2N  3 (N ≥ 2). Although there exist functional sites with spacers ≥ 15 bp, experiments suggest that their contribution to the overall binding affinity within a clustersite is negligible.
Figure 7. Comparison of Clustersite scores and Luciferase Activity. This graph compares the estimated relative binding affinity given by the clustersite score to the luciferase activity from four experiments for four different p53 clustersites. The four clustersites regulate the genes DDB2 (blue), CKM (red), IGFBP3 (green), and TP53I3 (cyan). In all four experiments the luciferase activity of truncated mutants of the respective p53 clustersite were compared to the luciferase activity of the full clustersite. In the case of the TP53I3 clustersite, four different mutants of varying lengths were measured for luciferase activity. All clustersite scores and activity measurements are normalized by the fullsite (two halfsites) measurement. The clustersite scores are attained by summing the estimated binding affinity of all viable fullsites in the clustersite that have an affinity above a lower bound and spacerlengths below an upper bound. The fullsite affinity lower bound and spacerlength upper bound were chosen to best match the experimental data. The best fit was attained by enforcing that spacerlengths not exceed 14 bp and affinity scores exceed 27.5.
Table 1. Normalized Experimental Affinity of Clustersites
These p53 clustersite scores are attained through a two step process. The first step uses the clustersite model which contains a generalized p53 halfsite PHMM and a backtransition that limits any spacer between two halfsites to no more than 14 bp (see part e of Figure 2). The dynamic programming Viterbi algorithm is used to find the highest scoring p53 halfsites in the sequence (that are separated by no more then 14 bp). The second step then parses the statepath generated from step 1 and generates viable p53 fullsites with any spacers removed, while conserving the property that the halfsites in the clustersite were not separated by more than 14 bp. Now we use the more flexible p53 singlesite model to score these viable fullsites using the Viterbi algorithm (see part d of Figure 2). We maintain a running sum of the logodds scores of the candidate fullsites that are above a certain threshold. The logodds score threshold and spacerlength limit (14 bp) are chosen so as to best fit the experimental data (see Figure 7).
Additionally, this p53 clustersite model follows statistical mechanics, in that the overall binding affinity for the complete RE is proportional to the probability of any p53 protein binding to any of the allowed motifs found in the clustersite. (See Methods for more details.)
Dynamic acceptance thresholds as a function of the distance from the TSS
An interesting finding from the analysis of our dataset of 160 functional p53 binding sites is that the low relative affinity scores from our model are significantly correlated with short distances from the Transcription Start Site (TSS). We find that low affinity sites exist only in a tight band around the TSS (see part a of Figure 8). Therefore a dynamic bindingaffinity acceptance threshold, dependent upon the putative site's distance from the TSS, can greatly reduce the false positive rate of our classifier. With a dynamic acceptance threshold, putative sites will require higher calculated binding affinities as their distance from the TSS increases in order to be accepted as potentially functional. For example, consider the linear dynamic acceptance threshold .00107·ΔX + 65.16 shown in Figure 8, with the additional restriction that the putative sites must be within 5,000 bp upstream and 1,000 bp downstream of the gene. Let the static acceptance threshold be all normalized scores above 70 with the same restriction that the putative sites must be within 5,000 bp upstream and 1,000 bp downstream of the gene. Even though the restricted dynamic threshold has a false negative rate of 22 out of 158 validated p53 sites (13.9%), and the restricted static threshold 32 out of 158 (20.3%), the restricted static threshold generates over 3.2 times as many positive hits when scoring all 39,288 isoforms of known genes in the human genome (hg18). Thus, the dynamic acceptance threshold has a lower known false negative rate and a considerably lower false positive rate. Different dynamic acceptance thresholds can be chosen to match desired levels of the known false negative rate and the genome hit rate (see part b of Figure 8). An important consideration when choosing an acceptance threshold is that a decrease in the threshold will in general produce an exponential increase in the number of positive hits.
Figure 8. Normalized affinity scores versus distances from the TSS. (Upper) This plot presents the normalized affinity scores returned from the p53 singlesite model versus the distance from the Transcription Start Site (TSS) for 158 experimentally validated p53binding sites. Low affinity sites exist in a tight band around the TSS (cyan vertical line). p53 activationsites are plotted in green, repressionsites in red, and both activation and repression in black. All sites ≥ 11 Kb from the TSS have relative affinity scores above the average of ≈ 78 (purple horizontal line). (Lower) This plot presents the estimated normalized affinity scores versus the positive distance (absolute value) from the TSS. Three dynamic acceptance thresholds are shown for scoring for putative p53 binding sites. The blue linear threshold corresponds to the formula .00107·ΔX + 65.16 and has a false negative rate of 18 out of 158 validated p53 sites (11.4%). The orange logarithmic threshold corresponds to the formula 9.6854·log(ΔX + 593.31)  15.308 and has a false negative rate of 5 out of 158 validated p53 sites (3.2%). Finally, the black squareroot threshold corresponds to the formula .23186·sqrt(ΔX + 7.5231) + 45.6 and has a false negative rate of 1 out of 158 validated p53 sites (0.63%). (ΔX = distance from TSS)
Conclusion
Profile Hidden Markov Models (PHMMs) can boost predictive power over weight matrices (PSSMs) when the binding motif is highly degenerative and tolerates insertions and/or deletions at various positions. The increase in predictive power for the p53binding motif can be seen in Figures 4 and 5. When the RE has a known repeated and/or palindromic motif, this prior knowledge can be used to correspond parameters in the model to exploit the redundancy in the information in the motif. We propose a novel "Corresponded BaumWelch" training algorithm that significantly boosts the predictive power of the p53RE model, as seen in Figures 4 and 5. When the motif is not known, all possible motifs for the given size can be sampled and crossvalidation techniques leveraged to infer the correct motif that maximizes predictive power. For example, Figure 5 reveals that the maximally predictive p53binding motif corresponds the four quartersites in a combinedpalindromic structure.
Our algorithms demonstrate the best predictive capability to date in classifying putative p53 binding sites. One algorithm uses a novel "Corresponded BaumWelch" training method that exploits the repeated, palindromic structure of the p53 motif to train for allowed insertions and deletions relative to the consensus. The second algorithm properly models the relative increase in binding affinity for p53 clustersites (REs with ≥ 3 adjacent halfsites) by using a two step process that scores all viable fullsites in the clustersite while restricting the spacerlength to 14 bp. This new clustersite algorithm best matches the experimental data (see Figure 7).
Functional lowaffinity p53sites only exist near the TSS. Therefore the binding affinity threshold for accepting a putative site should be dependent on the putative site's distance from the TSS. By this method, putative sites with relatively low calculated binding affinities that are near the TSS may be accepted, while those sites with equal scores but more distant from the TSS will be rejected. A dynamic threshold, as a function of the distance from the TSS, can greatly reduce the false positive rate when searching for putative p53sites in genes.
Methods
The Corresponded BaumWelch algorithm
In order to exploit the redundancy of information in a homodimer or homotetramer binding motif, we wish to share information between corresponding positions. The level of sharing of information for any set of corresponding positions is given by a correspondence factor c such that 0 ≤ c ≤ 1. At the end of each round of the iterative BaumWelch algorithm we calculate the average values of each of the newly updated emission probabilities and transition probabilities for all k and l in the set of corresponding positions, represented as and respectively. Each of these average values represents the expected probability if the corresponding positions are fully tied (c = 1), and are referred to as the "corresponding average". Then we update the new emission and transition probabilities within the set of corresponding positions, using the current correspondence factor and corresponding average, according to:
If we wish to train for the optimum correspondence factor, then we calculate a new c' for each emission and transition probability at each position in the set of corresponding positions:
Now, we can calculate a new correspondence factor c' by averaging over sets of the and values. The one optimum correspondence factor for the whole motif or separate correspondence factors for sets of corresponding positions are obtained by averaging over different sets:
The Corresponded BaumWelch algorithm will converge at (local) optimum emission and transition probabilities and correspondence factors that maximize the likelihood of observing the training set with possible pseudocounts. Please see the Additional file 1 for further details.
The p53 clustersite algorithm
The p53 clustersite algorithm is a two step process designed to sum the estimated relative binding affinities of all viable fullsites within a clustersite. The first step uses the clustersite model that contains a generalized p53 halfsite PHMM and a backtransition through a nocost FEM14 module (see part e of Figure 2). The nocost Finite Emission Module (FEM) of length 14 can match any sequence of length ≤ 14 bp with no contribution to the overall score. We score the entire putative clustersite using the p53 clustersite model and the Viterbi algorithm to find the bestsupported path through the clustersite. This path provides the strongest affinity halfsites that are not separated by more than 14 bp. If we use the notation "14" for any spacer sequence of length 0 to 14 and H for a halfsite sequence, then we can represent the clustersite sequence path as:
Step 2 now parses the clustersite sequence path and generates a list of all viable fullsites, which are concatenations of any two halfsites such that they are not separated by more than 14 bp:
Now we use the more flexible (and more accurate) singlesite model with the Viterbi algorithm to estimate the relative binding affinity of all the viable fullsites in the clustersite. The clustersite affinity score is the sum of all viable fullsite scores that exceed a certain threshold. If F denotes a viable fullsite then:
The spacerlength upper bound and the affinityscore lower bound were fit to best match the experimental results. In the case for p53binding sites, the best fit is a spacerlength of no more than 14 bp and a logodds score of at least 27.5 (see Figure 7).
The p53HMM implementation
The p53HMM algorithm is implemented in Java and is available online at http://tools.csb.ias.edu webcite. The implementation makes extensive use of the BioJava Toolkit [24].
Authors' contributions
TR participated in the design of the algorithms, wrote the code, performed the computational analysis, and drafted the manuscript. XY performed all experiments. ES participated in the design of the algorithms and helped to draft the manuscript. AL conceived of the study and helped to draft the manuscript. All authors read, edited, and approved the final version of the paper.
Acknowledgements
We thank Michael Krasnitz, Amar Drawid, Anirvan Sengupta, Sean Eddy, Jiri Vanicek, and Raúl Rabadán for helpful discussions.
References

Levine AJ: p53, the cellular gatekeeper for growth and division.
Cell 1997, 88(3):323331. PubMed Abstract  Publisher Full Text

Funk WD, Pak DT, Karas RH, Wright WE, Shay JW: A transcriptionally active DNAbinding site for human p53 protein complexes.
Mol Cell Biol 1992, 12(6):28662871. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

el Deiry WS, Kern SE, Pietenpol JA, Kinzler KW, Vogelstein B: Definition of a consensus binding site for p53.
Nat Genet 1992, 1:4549. PubMed Abstract  Publisher Full Text

Riley T, Sontag E, Chen P, Levine A: Transcriptional control of human p53regulated genes.
Nat Rev Mol Cell Biol 2008, 9(5):402412. PubMed Abstract  Publisher Full Text

Krogh A, Brown M, Mian IS, Sjölander K, Haussler D: Hidden Markov models in computational biology. Applications to protein modeling.
J Mol Biol 1994, 235(5):15011531. PubMed Abstract  Publisher Full Text

Eddy SR: Profile hidden Markov models.
Bioinformatics 1998, 14(9):755763. PubMed Abstract  Publisher Full Text

Stormo G, Fields D: Specificity, free energy and information content in proteinDNA interactions.
Trends in Biochemical Sciences 1998, 23(5):109113. PubMed Abstract  Publisher Full Text

Durbin R, Eddy S, Krogh A, Mitchison G: Biological sequence analysis. 1st edition. Cambridge University Press; 1998.

Djordjevic M, Sengupta AM, Shraiman BI: A Biophysical Approach to Transcription Factor Binding Site Discovery.
Genome Res 2003, 13(11):23812390. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lee J, Kim J, Kim J: Datadriven design of hmm topology for online handwriting recognition. [http://citeseer.ist.psu.edu/lee00datadriven.html] webcite
2000.

Ma B, Pan Y, Zheng J, Levine AJ, Nussinov R: Sequence analysis of p53 responseelements suggests multiple binding modes of the p53 tetramer to DNA targets.
Nucleic Acids Res 2007, 35(9):29863001. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Marinescu VD, Kohane IS, Riva A: MAPPER: a search engine for the computational identification of putative transcription factor binding sites in multiple genomes.
BMC Bioinformatics 2005, 6:79. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yu X, Riley T, Levine AJ: The regulation of the endosomal compartment by p53 the tumor suppressor gene.
FEBS Journal 2009, 276(8):22012212. PubMed Abstract  Publisher Full Text

Hu W, Feng Z, Teresky AK, Levine AJ: p53 regulates maternal reproduction through LIF.
Nature 2007, 450(7170):721724. PubMed Abstract  Publisher Full Text

Hoh J, Jin S, Parrado T, Edington J, Levine AJ, Ott J: The p53MH algorithm and its application in detecting p53responsive genes.
Proc Natl Acad Sci USA 2002, 99(13):84678472. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hughey R, Krogh A: Hidden Markov models for sequence analysis: extension and analysis of the basic method.
Comput Appl Biosci 1996, 12(2):95107. PubMed Abstract

Barrett C, Hughey R, Karplus K: Scoring hidden Markov models.
Comput Appl Biosci 1997, 13(2):191199. PubMed Abstract  Publisher Full Text

Barash Y, Elidan G, Friedman N, Kaplan T: Modeling dependencies in proteinDNA binding sites. In RECOMB '03: Proceedings of the seventh annual international conference on Research in computational molecular biology. New York, NY, USA: ACM Press; 2003:2837.

Zhou Q, Liu JS: Modeling withinmotif dependence for transcription factor binding site predictions.
Bioinformatics 2004, 20(6):909916. PubMed Abstract  Publisher Full Text

Tan T, Chu G: p53 Binds and activates the xeroderma pigmentosum DDB2 gene in humans but not mice.
Mol Cell Biol 2002, 22(10):32473254. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Contente A, Dittmer A, Koch MC, Roth J, Dobbelstein M: A polymorphic microsatellite that mediates induction of PIG3 by p53.
Nat Genet 2002, 30(3):315320. PubMed Abstract  Publisher Full Text

Bourdon JC, DeguinChambon V, Lelong JC, Dessen P, May P, Debuire B, May E: Further characterisation of the p53 responsive elementidentification of new candidate genes for transactivation by p53.
Oncogene 1997, 14:8594. PubMed Abstract  Publisher Full Text

Kern SE, Pietenpol JA, Thiagalingam S, Seymour A, Kinzler KW, Vogelstein B: Oncogenic forms of p53 inhibit p53regulated gene expression.
Science 1992, 256(5058):827830. PubMed Abstract  Publisher Full Text

Holland RCG, Down TA, Pocock M, Prlić A, Huen D, James K, Foisy S, Dräger A, Yates A, Heuer M, Schreiber MJ: BioJava: an opensource framework for bioinformatics.
Bioinformatics 2008, 24(18):20962097. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Schneider TD, Stephens RM: Sequence logos: a new way to display consensus sequences.
Nucleic Acids Res 1990, 18(20):60976100. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

SchusterBöckler B, Schultz J, Rahmann S: HMM Logos for visualization of protein families.
BMC Bioinformatics 2004, 5:7. PubMed Abstract  Publisher Full Text  PubMed Central Full Text