- Research article
- Open access
- Published:

# The p53HMM algorithm: using profile hidden markov models to detect p53-responsive genes

*BMC Bioinformatics*
**volume 10**, Article number: 111 (2009)

## 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 single-sites and cluster-sites. These models incorporate a novel "Corresponded Baum-Welch" training algorithm that provides increased predictive power by exploiting the redundancy of information found in the repeated, palindromic p53-binding 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 p53-binding sites increases predictive accuracy by reducing the false positive rate.

### Results

Training a Profile Hidden Markov Model with corresponding positions matching a combined-palindromic p53-binding motif creates the best p53-RE predictive model. The p53HMM algorithm is available on-line: http://tools.csb.ias.edu

### 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 cancer-inducing, DNA-damaging 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 half-site RRRCWWGYYY, followed by a variable length spacer, then followed (almost always) by a second half-site 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 quarter-site RRRCW as → and its reverse-complement WGYYY as ←, the first discovered p53 consensus sequence can be graphically represented by → ← *spacer* → ←. This configuration of the four quarter-sites is often referred to as the head-to-head (HH) orientation, and represents the vast majority of experimentally-validated p53 binding sites to date.

### The degeneracy of the p53-RE

In the influential paper "Definition of a Consensus binding Site for p53", by El-Deiry et al., 7 of the 20 DNA target sites (35%) used to form the head-to-head (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 experimentally-validated 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 mis-score 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 non-critical nucleotides (including tolerated insertions in the motif) are not conserved. The match-state 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 fine-tune 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).

### 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 log-odds 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 TF-protein binding to that candidate site [see Additional file 1 for details] [7–9]. The log-odds 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 binding-site 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 log-odds score for the sequence [8]. In the case of transcription factor binding sites, the log-odds 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 protein-DNA complex, and that all other possible binding locations in the response element can be ignored.

### 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 over-fitting. Again, under ideal conditions the log-odds score (*log-likelihood 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 *Baum-Welch* algorithm is the most commonly used iterative Expectation Maximization (EM) method to train the parameters of the model. The Baum-Welch 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 p53-motif 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 search-space 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 protein-DNA 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 combined-palindromic motif that completely corresponds the four quarter-sites is the *discovered* motif, since it is the best classifier (see Figure 4).

### The Corresponded Baum-Welch algorithm

In order to include the prior knowledge of the structural motif (or in an attempt to discover it), a novel "Corresponded Baum-Welch" algorithm is proposed to enforce or learn the optimal correspondence between expectations of parameters for corresponding positions after each iteration of the Baum-Welch 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 reverse-complement 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 bp-tetrameric 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 T-coupled motif (*1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10*), the (reverse-complement) palindromic or H-coupled motif (*1 2 3 4 5 6 7 8 9 10* ), the independently (reverse-complement) palindromic or Q-coupled motif (*1 2 3 4 5* *6 7* *8 9 10* , the repeated, fully-palindromic or combined-palindromic motif (*1 2 3 4 5* *1 2 3 4 5* ), and the completely un-tied motif with no correspondence between any positions (see Figure 3) [11]. We perform 1000 iterations of ten-fold random-split 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, threshold-independent (non-parametric) 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 Insert-State 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 over-all average. The more non-fixed 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 insertion-state 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 insertion-state emission distributions at every position to the amino-acid average for the whole sequence would be very inappropriate in this example.

### The final results

The combined-palindromic 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 combined-palindromic motif leverages the correlations found in all of them. Furthermore, it can be seen that training the insert-state 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 insert-state emission training to the base PHMM (with no motif-corresponded positions) has an AUC improvement of .923 - .919 = .004, but with motif training has one of .937 - .929 = .008. Adding motif training (motif-corresponded positions) to the PHMM when not insert-state emission training has an AUC improvement of .929 - .919 = .010, but with insert-state emission training has one of .937 - .923 = .014. Therefore the improvements are not additive. There is "positive synergy" when performing both motif training and insert-state emission training together that further boosts the predictive accuracy of the model. This observation confirms our hypothesis that training insert-state emissions can significantly boost the accuracy of the model after corresponding positions in the PHMM according to a binding-site motif.

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 combined-palindromic 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 quarter-sites in a combined-palindromic 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 combined-palindromic structural motif, and training the insert-state emissions (see Figure 6).

### 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 p53-regulated 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 half-site relative to the recently published putative site in intron 1 [14]. The p53HMM algorithm predicted the site GGACATGTCGGGACA-GCTC, which matches the consensus RRRCWWGYYYRRRCWWGYYY perfectly except for the low-conserved position 10 and the gap ("-", deletion) at position 16. A PSSM approach predicted the shifted site GGACATGTCGggacagCTCCCAGCTC, which is the best "gap-less" 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 well-conserved 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 arbitrary-length spacers by inserting a no-cost *Free Insertion Module* (FIM) between the two half-sites of the single-site PHMM [16, 17]. Similarly, we can ignore spacers with lengths between 1 and *N* base pairs by inserting a no-cost *Finite Emission Module* (FEM-N) between the two half-sites (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 over-fit 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 reverse-complement received a score of 0, and vice-versa. 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 over-fitting 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 protein-DNA 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 higher-order 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 cluster-site algorithm

Binding affinity measurements have been obtained for certain p53 cluster-sites of different lengths by mutating or truncating known p53 cluster-sites in the genes: DDB2, TP53i3, CKM, IGFBP3, and RGC (see Table 1 and Figure 7) [20–23]. Based on the relative binding affinities of these p53 cluster-sites, we propose a new p53 cluster-site algorithm that utilizes the trained PHMM to calculate and sum up the relative estimated binding-affinities, above a certain threshold, of all viable full-sites 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 half-sites in the cluster-site and the length of spacers between them. For example, for p53 cluster-sites with 2, 3, 4, 5, or 6 adjacent p53 half-sites, the number of possible full-sites with spacer-lengths = 14 bp would be 1, 3, 5, 7, and 9, respectively. Let N be the number of half-sites in the cluster-site, then the number of full-sites (to calculate binding affinities for and sum up) is given by the expression 2*N* - 3 (*N* ≥ 2). Although there exist functional sites with spacers ≥ 15 bp, experiments suggest that their contribution to the overall binding affinity within a cluster-site is negligible.

These p53 cluster-site scores are attained through a two step process. The first step uses the cluster-site model which contains a generalized p53 half-site PHMM and a back-transition that limits any spacer between two half-sites 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 half-sites in the sequence (that are separated by no more then 14 bp). The second step then parses the state-path generated from step 1 and generates viable p53 full-sites with any spacers removed, while conserving the property that the half-sites in the cluster-site were not separated by more than 14 bp. Now we use the more flexible p53 single-site model to score these viable full-sites using the *Viterbi* algorithm (see part d of Figure 2). We maintain a running sum of the log-odds scores of the candidate full-sites that are above a certain threshold. The log-odds score threshold and spacer-length limit (14 bp) are chosen so as to best fit the experimental data (see Figure 7).

Additionally, this p53 cluster-site 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 cluster-site. (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 binding-affinity 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.

## 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 p53-binding 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 Baum-Welch" training algorithm that significantly boosts the predictive power of the p53-RE 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 cross-validation techniques leveraged to infer the correct motif that maximizes predictive power. For example, Figure 5 reveals that the maximally predictive p53-binding motif corresponds the four quarter-sites in a combined-palindromic structure.

Our algorithms demonstrate the best predictive capability to date in classifying putative p53 binding sites. One algorithm uses a novel "Corresponded Baum-Welch" 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 cluster-sites (REs with ≥ 3 adjacent half-sites) by using a two step process that scores all viable full-sites in the cluster-site while restricting the spacer-length to 14 bp. This new cluster-site algorithm best matches the experimental data (see Figure 7).

Functional low-affinity p53-sites 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 p53-sites in genes.

## Methods

### The Corresponded Baum-Welch 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 Baum-Welch 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 Baum-Welch algorithm will converge at (local) optimum emission and transition probabilities and correspondence factors that maximize the likelihood of observing the training set with possible pseudo-counts. Please see the Additional file 1 for further details.

### The p53 cluster-site algorithm

The p53 cluster-site algorithm is a two step process designed to sum the estimated relative binding affinities of all viable full-sites within a cluster-site. The first step uses the cluster-site model that contains a generalized p53 half-site PHMM and a back-transition through a no-cost FEM-14 module (see part e of Figure 2). The no-cost Finite Emission Module (FEM) of length 14 can match any sequence of length ≤ 14 bp with no contribution to the over-all score. We score the entire putative cluster-site using the p53 cluster-site model and the *Viterbi* algorithm to find the best-supported path through the cluster-site. This path provides the strongest affinity half-sites 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 half-site sequence, then we can represent the cluster-site sequence path as:

Step 2 now parses the cluster-site sequence path and generates a list of all viable full-sites, which are concatenations of any two half-sites such that they are not separated by more than 14 bp:

Now we use the more flexible (and more accurate) single-site model with the Viterbi algorithm to estimate the relative binding affinity of all the viable full-sites in the cluster-site. The cluster-site affinity score is the sum of all viable full-site scores that exceed a certain threshold. If *F* denotes a viable full-site then:

The spacer-length upper bound and the affinity-score lower bound were fit to best match the experimental results. In the case for p53-binding sites, the best fit is a spacer-length of no more than 14 bp and a log-odds score of at least 27.5 (see Figure 7).

### The p53HMM implementation

The p53HMM algorithm is implemented in Java and is available on-line at http://tools.csb.ias.edu. The implementation makes extensive use of the BioJava Toolkit [24].

## References

Levine AJ: p53, the cellular gatekeeper for growth and division.

*Cell*1997, 88(3):323–331.Funk WD, Pak DT, Karas RH, Wright WE, Shay JW: A transcriptionally active DNA-binding site for human p53 protein complexes.

*Mol Cell Biol*1992, 12(6):2866–2871.el Deiry WS, Kern SE, Pietenpol JA, Kinzler KW, Vogelstein B: Definition of a consensus binding site for p53.

*Nat Genet*1992, 1: 45–49.Riley T, Sontag E, Chen P, Levine A: Transcriptional control of human p53-regulated genes.

*Nat Rev Mol Cell Biol*2008, 9(5):402–412.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):1501–1531.Eddy SR: Profile hidden Markov models.

*Bioinformatics*1998, 14(9):755–763.Stormo G, Fields D: Specificity, free energy and information content in protein-DNA interactions.

*Trends in Biochemical Sciences*1998, 23(5):109–113.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):2381–2390.Lee J, Kim J, Kim J: Data-driven design of hmm topology for on-line handwriting recognition.2000. [http://citeseer.ist.psu.edu/lee00datadriven.html]

Ma B, Pan Y, Zheng J, Levine AJ, Nussinov R: Sequence analysis of p53 response-elements suggests multiple binding modes of the p53 tetramer to DNA targets.

*Nucleic Acids Res*2007, 35(9):2986–3001.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.Yu X, Riley T, Levine AJ: The regulation of the endosomal compartment by p53 the tumor suppressor gene.

*FEBS Journal*2009, 276(8):2201–2212.Hu W, Feng Z, Teresky AK, Levine AJ: p53 regulates maternal reproduction through LIF.

*Nature*2007, 450(7170):721–724.Hoh J, Jin S, Parrado T, Edington J, Levine AJ, Ott J: The p53MH algorithm and its application in detecting p53-responsive genes.

*Proc Natl Acad Sci USA*2002, 99(13):8467–8472.Hughey R, Krogh A: Hidden Markov models for sequence analysis: extension and analysis of the basic method.

*Comput Appl Biosci*1996, 12(2):95–107.Barrett C, Hughey R, Karplus K: Scoring hidden Markov models.

*Comput Appl Biosci*1997, 13(2):191–199.Barash Y, Elidan G, Friedman N, Kaplan T: Modeling dependencies in protein-DNA 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:28–37.Zhou Q, Liu JS: Modeling within-motif dependence for transcription factor binding site predictions.

*Bioinformatics*2004, 20(6):909–916.Tan T, Chu G: p53 Binds and activates the xeroderma pigmentosum DDB2 gene in humans but not mice.

*Mol Cell Biol*2002, 22(10):3247–3254.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):315–320.Bourdon JC, Deguin-Chambon V, Lelong JC, Dessen P, May P, Debuire B, May E: Further characterisation of the p53 responsive element-identification of new candidate genes for trans-activation by p53.

*Oncogene*1997, 14: 85–94.Kern SE, Pietenpol JA, Thiagalingam S, Seymour A, Kinzler KW, Vogelstein B: Oncogenic forms of p53 inhibit p53-regulated gene expression.

*Science*1992, 256(5058):827–830.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 open-source framework for bioinformatics.

*Bioinformatics*2008, 24(18):2096–2097.Schneider TD, Stephens RM: Sequence logos: a new way to display consensus sequences.

*Nucleic Acids Res*1990, 18(20):6097–6100.Schuster-Böckler B, Schultz J, Rahmann S: HMM Logos for visualization of protein families.

*BMC Bioinformatics*2004, 5: 7.

## Acknowledgements

We thank Michael Krasnitz, Amar Drawid, Anirvan Sengupta, Sean Eddy, Jiri Vanicek, and Raúl Rabadán for helpful discussions.

## Author information

### Authors and Affiliations

### Corresponding author

## Additional information

### 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.

## Electronic supplementary material

### 12859_2008_2841_MOESM1_ESM.pdf

Additional File 1: **Supplementary Material**. The Supplementary Material contains more theory behind modeling TF-binding sites with PHMMs, details of the Corresponded Baum-Welch algorithm, and a proof that the PHMM log-odds score of a TF-binding site estimates its relative binding affinity given certain assumptions. (PDF 201 KB)

## Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

## Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## About this article

### Cite this article

Riley, T., Yu, X., Sontag, E. *et al.* The p53HMM algorithm: using profile hidden markov models to detect p53-responsive genes.
*BMC Bioinformatics* **10**, 111 (2009). https://doi.org/10.1186/1471-2105-10-111

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/1471-2105-10-111