BACKGROUND
identification of short, degenerate dna sequences  that bind to a transcription factor is a difficult problem  <cit> .
the particular sub-problems of identification of sites corresponding to self-overlapping motifs, determination of a threshold for this purpose, biophysical interpretation of statistical quantities used in probabilistic site identification methods and estimation of occupancy of these sites by a transcription factor have not been addressed in detail before.
these problems form the focus of this study.
several transcription factors bind to self-overlapping dna motifs.
when a motif overlaps with itself with a particular shift, the transcription factor can bind in more than one sequence windows, i.e.
to overlapping sites.
examples include drosophila developmental transcription factor hunchback  <cit> , worm pha- <dig>  <cit> , human sp- <dig>  c/ebpalpha, yeast adr <dig>  mig <dig>  chicken cdx- <dig>  arabidopsis agamous, etc.
<cit> .
furthermore, when binding by the transcription factor in either orientation is permissible, the corresponding dna site and its reverse complement can be considered as two different types of sites that are overlapping, as long as the motif is not exactly palindromic.
problems similar to the identification of overlapping sites also arise in the context of prediction of nucleosome positioning  <cit> .
self-overlapping motif of nf-κb family of factors
the nf-κb family of transcription factors is a prominent example of transcription factors with self-overlapping motifs.
these factors are key mediators of the cellular response to infection, injury, inflammation, or other stress conditions that lead to rapid alterations in cellular gene expression  <cit> .
the vertebrate members include the c-rel, rela, relb, nfkb1/p <dig> and nfkb2/p <dig> proteins that form homo- or heterodimers with one another and bind to similar sites.
these nf-κb binding sites , with the consensus being gggrnnyycc  <cit> , often overlap because a κb site contains multiple g's at the 5' end and multiple c's at the 3' end with a high probability .
for a good κb site, when the sequence window is shifted by one base in either the 5' or the 3' direction the resulting sequence is also a putative κb site.
moreover, we can have additional possibility of self-overlap because the reverse complement of a κb site is often a κb site, allowing functional binding in the opposite direction.
many computational methods do not take into account these alternative binding modes while scoring a candidate regulatory site.
ignoring the self-overlapping nature of the motif can also complicate the interpretation of in vitro experiments.
for example, a 3' padding sequence, starting with nucleotide c, in a gel shift experiment can form a spurious strong κb site and thus confer higher binding affinity to the experimental sequence even when the test sequence, in the context of the native promoter, makes for a very weak κb site.
a computational method that naturally considers multiple binding windows and help in the design of padding sequences is useful for avoiding this problem.
the commonly used site identification methods often assume equal probability of a site everywhere in a particular window around the transcription start site  and thus fail to adjust according to the observed distribution of sites within the gene structure.
whereas proximal promoters, up to  <dig> bp upstream of the tss, contain several sites, further upstream distal promoters contain fewer sites  <cit> .
while identifying sites, a method should adjust its parameters to reflect this change in the density of sites with respect to the distance from the tss.
a computational method that links to the biophysical process of transcriptional regulation can offer more biological insight than most of the current site identification methods  <cit> .
such a method can interpret statistical quantities as biophysical variables like binding energy and transcription factor concentration.
it can also estimate how often a transcription factor is bound to the gene promoter, i.e.
its occupancy probability.
interestingly, the use of occupancy probability as a discriminant function during site identification has a number of advantages .
one advantage of this approach is that it goes beyond just site identification and provides a measure of the importance of a site in modulating gene regulation.
for example, for an activator, highly occupied sites may result in high gene expression and thus may be biologically more significant.
most methods that identify transcription factor binding sites or regulatory modules, however, leave unexplored the occupancy of the transcription factor on the promoter and its impact on gene regulation.
the lack of biophysical interpretation prevents these methods from using a natural threshold for classifying sequences into sites.
they declare that a sequence is a site if its score is above an arbitrary threshold or is statistically significant compared to the score of the random background.
they often have a high threshold resulting in several false negatives or a low threshold resulting in several false positives.
thus, identification of weak sites is difficult.
on the contrary, saturation effects in the physical occupancy provides a natural threshold and this fact has been utilized to develop class support vector machines like qpmeme  <cit> .
alternatively, hidden markov models , when interpreted as biophysical models, overcome this problem by learning the threshold based on the transition probability to the motif state, namely the hidden state associated with a sequence pattern which is distinct from the background.
hmm for binding site identification
hmms, now popular in bioinformatics for more than a decade  <cit> , have been used in two different ways for motif identification.
for identification of one or more occurrences of non-overlapping sites: 'profile hmms'  <cit> , originally designed to model protein domains, have more recently been used to identify binding sites of transcription factors, for example, of camp receptor protein in cyanobacterium anabaena  <cit> , liver x receptor  <cit>  and creb  <cit> .
a profile hmm library was built using transfac sequences to classify transcription factors  <cit> .
in a profile hmm, each position within a motif has three states.
a match state is associated with a nucleotide being present at that position and has corresponding emission probabilities.
a deletion state corresponds to absence of any nucleotide at that position.
an insertion state allows for insertion of nucleotides between the current position and the next position within the motif, and has its own emission probabilities.
for identification of cis-regulatory modules  that contain multiple sites of different types: this is usually performed using 'motif hmms'  <cit> .
in a motif hmm, the entire motif, i.e.
all positions within the motif, is represented by one state.
different states correspond to different motif types .
phylogenetic conservation has been incorporated in motif hmms to reduce false positives  <cit> .
both types of hmms, however, do not generally train their parameters properly.
a profile hmm has a complicated architecture and requires a large number of parameters as a consequence.
because the number of known sites is generally small, training of a profile hmm using known sites in their native promoters  is usually not possible.
therefore, a profile hmm is generally trained using pre-aligned sites.
because the whole promoters containing the training sites are not used, transition probability to the motif z from background is not trained, and the relationship between transition and emission probabilities is likewise not captured.
motif hmms, on the other hand, are more focused on identifying motifs of multiple types.
but while they attempt to estimate the transition probability to the motif z, they generally use only the promoters containing the motifs for training z and thus overestimate z. moreover, they usually train emission probabilities of motifs separately using pre-aligned training sites, thus ignoring the effect of z on emission probabilities.
many of these hmms use the likelihood or the viterbi algorithm for scoring , and thus end up using an arbitrary classification threshold.
furthermore, they leave their relationship with biophysical models rather obscure and thus fail to calculate the occupancy probability.
these shortcomings are in addition to their basic failure to explicitly consider the effect of multiple overlapping sites in the viterbi method.
methods specifically developed for κb site identification
in an original approach to specifically identify κb sites, udalova et al.
developed a principal coordinate model  <cit> .
they determined relative binding affinities of κb sites using experimental quantitative binding data.
they selected a subset of the  <dig> possible variants of the fully palindromic nf-κb binding consensus sequence ggrrnnyycc such that no variant differed from the selected sequences or their reverse complements by more than one nucleotide.
they mapped these sequences to a euclidean space and used the largest principal components for least-square linear regression of the logarithm of binding affinity in a gel shift assay or microarray.
this model automatically incorporated effects of interactions between base pair positions in the binding site, its predictions were highly correlated with experimental binding data and it identified site positions responsible for differential binding of nf-κb family members consistent with crystallographic studies.
however, the model's disadvantages are that it  requires experimental quantitative binding data of all selected sequences and  only includes variants of the consensus sequence.
because several known κb sites do not fit the consensus sequence  <cit> , inclusion of all possible 10-mer variants for this model will require binding experiments with a large number of sequences, making this model infeasible.
furthermore, this model also suffers from the limitations of site identifications described above.
alternative interpretation of hmm as a physical binding model
even though in the context of site identification an hmm is usually interpreted as a sequence generative model, we focus on a less familiar interpretation as a physical binding model of a transcription factor on dna.
this interpretation leads us to transform the statistical hmm model into a biophysical one.
we can then determine the occupancy probability of a transcription factor on a dna sequence, and think of the prior probability of the motif  as a measure of transcription factor concentration and the weight matrix  as a measure of binding energies.
more importantly, the biophysical model offers a principled threshold for classifying sequences into sites.
the statistical mechanics model was utilized in  <cit>  while the equivalent hmm representation is used in  <cit> .
although we do not know of a paper explicitly describing the equivalence of these two approaches, the work on nucleosome positioning by segal et al.
<cit>  utilizes this relationship.
an hmm is commonly used as a generative model of a spatial or time series sequence in the machine learning field  <cit> .
in our context, the hmm consisting of two kinds of states – the background state and motif states – generates a sequence from 5' to 3'.
at any position in the sequence, the hmm  determines the probabilities of the motif and background states at the previous position,  calculates the probability of either state at the current position using the transition probabilities and  generates the new nucleotide based upon the hidden states' emission probabilities.
we could also consider the hmm as a physical binding model that estimates the occupancy probability of a transcription factor on a particular position of a dna sequence, i.e.
how often the transcription factor is bound to that position of the dna sequence.
the relation between the two interpretations of an hmm, are as follows.
in our particular hmm, the start positions of a motif state, with a particular orientation, specifies the hidden states.
one could ask define the conditional probability distribution of the hidden states, or, equivalently these positions, is given the parameters and the observed sequence.
it turns out that this conditional distribution is exactly the same as that of the edge positions of bound transcription factors on dna in the statistical mechanical model.
the occupancy probability at a position is the probability of finding a motif state at that position in the hmm language.
if we ignore the possibility of binding in multiple overlapping windows, then, the probability at the jth position of sequence s being the beginning of a motif is , where wj is the weight matrix score of the motif starting at the jth position of the sequence .
we claim that this quantity can be mapped to the occupancy probability, as we will see below.
thus, the two factors determining occupancy probability are  the transition probability to the motif z and  the measure of distinctness of the emission probabilities of the motif from that of the background .
both these factors need to be high for the transcription factor to be bound to a particular position in a dna sequence with high probability.
for example, even if the weight matrix score is high, occupancy probability cannot reach one if z is really small.
a site identification method based only on a weight matrix has no way of dealing with this interplay with z.
the hmm training techniques offer two advantages over the calculations made using first principles.
first, the hmm baum-welch procedure trains the transition and emission probabilities.
this produces optimized values of z and the weight matrix, which are essential for an accurate estimation of occupancy probability.
secondly, as mentioned above, the hmm training procedure also calculates an intermediate variable called the gamma  variable.
the occupancy probability at a particular position is given simply by the gamma variable of the motif state at that position: , where  and  are the gamma values of the motif and background states at that position, respectively .
thus, the calculation of occupancy probability in a new sequence requires a simple extension of the scoring procedure in which the gamma variable is computed just like during the training procedure.
in a biophysical model, based on thermodynamics principles, occupancy probability can be written as , where  is the concentration of a free transcription factor at equilibrium, ej is the binding energy at position j,  where kb is the boltzmann's constant, and t is the absolute temperature  <cit> .
a comparison of this equation with the equation of probability obtained from the hmm gives us wj = - βej and z = .
thus, the weight matrix represents binding energy  <cit> .
in addition, the transition probability to the motif z corresponds to the free transcription factor concentration.
as the transcription factor concentration increases, the transition probability to the motif increases and we expect higher occupancy by the transcription factor on the dna.
one main focus of this study is to determine a good threshold for the classification of sequences into binding sites.
use of an arbitrary threshold while using a purely statistical quantity as a discriminant function is prone to a high number of false positives or false negatives.
a good discriminant function has a physical interpretation, it relates to the biological significance of a sequence and it offers a natural threshold for classification.
moreover, the classification method needs to be able define the above threshold based on training sequences.
we use occupancy probability as a discriminant function to classify sequences into sites because it has the following benefits.
first, unlike a purely statistical entity, it has a straightforward physical interpretation.
a transcription factor's occupancy on the promoter determines gene expression.
second, highly occupied binding sites may be physiologically more significant.
thus, occupancy probability not only helps in classifying sequences as binding sites but also offers insight into their influence on gene expression.
third, occupancy probability has a "natural" threshold at  <dig>  due to its fermi-dirac distribution with respect to binding energy  <cit> .
sequences with binding energy less than the chemical potential have occupancy probability greater than  <dig>  and hence can be classified as sites.
because of our focus on occupancy, we call this model ohmm .
the greatest benefit of the transformation of an hmm to a biophysical model is that it enables the hmm to learn the threshold for classification of sequences into sites in a principled way.
we have seen above that the occupancy probability depends upon the transition probability to the motif z and the weight matrix, both of which are trained by an hmm.
thus, when an hmm uses occupancy probability as a discriminant function, it learns the natural threshold based on the training sequences.
because of this accurate estimation of the threshold, an hmm is expected to identify weak sites much more accurately with a fewer false positives.
in summary, we describe an hmm, applied to nf-κb, that identifies overlapping sites, learns the threshold in a principled manner by training emission probabilities using known sites in their native promoters and training transition probabilities using promoters in an entire genome, estimates the transcription factor's occupancy probability, and adjusts parameters reflecting the distance from the tss while scoring.
the structure of hidden states in ohmm is indicated in figure  <dig>  we demonstrate that the occupancy probability predicted by ohmm correlates well with the binding affinity in gel shift experiments.
high evolutionary conservation scores and enrichment in experimentally regulated genes suggest that the predicted sites might be functional.
RESULTS
hmm training
training ohmm consists of estimating two sets of parameters:  the emission probabilities of each motif state  and the background, and  the transition probabilities that depend upon the transition probability to the motif .
while an hmm needs to be trained using site-rich sequences to learn the motif profile, training on random sequences is required to learn z reflecting the site density in the promoters of all genes in the human genome.
we therefore estimated hmm parameters in two steps.
we first trained all hmm parameters using site-rich sequences and thus learned the emission probabilities.
we then fixed the emission probabilities and trained the transition probabilities using randomly selected human promoters to estimate z separately for the tss-800:tss  and tss:tss+ <dig>  regions, where tss is the transcription start site.
the reason for separate z estimation in these two types of regions is that κb site density is different upstream and downstream of the tss.
training the emission probabilities
in the first training step, we used two types of site-rich sequences of different lengths and various initial z's to train all hmm parameters and determine the emission probabilities to be used in further analysis.
trained motif profiles of "tss-n hmms," i.e.
hmms trained on promoters consisting of n bases upstream of the tss of human genes known to contain a κb site, appear similar to the background emission probabilities regardless of the promoter length and initial z. on the other hand, trained motif profiles of "surround-50" or "surround- <dig> hmms," i.e.
hmms trained on  <dig> or  <dig> bp sequences each consisting of a known κb site and surrounding region, with a reasonable initial z , are distinct from the background .
they are also distinct from the initial motif profile, as their symmetrical kullback-leibler  divergences  with respect to the initial motif profile are high .
the trained motif profiles are slightly weaker than the initial motif profile, i.e.
more similar to the background.
in a surround- <dig> or surround- <dig> hmm, any initial z between  <dig>  and  <dig>  results in the same trained motif profile, indicating that perhaps a local optimum is reached.
trained motif profiles of surround- <dig> hmms, however, appear more and more like the background as the initial z increases above  <dig> .
trained motif profiles of surround- <dig> hmms appear similar to the background regardless of the initial z. we used the trained motif profile of the surround- <dig> hmm for further analysis .
estimation of the transition probability to the motif
in the second training step, we estimated z by training the transition probabilities of the above surround- <dig> hmm using random sequences while keeping its emission probabilities constant.
we estimated z separately for the upstream  <dig> bp and downstream  <dig> bp regions.
in each case, we determined human promoters with known κb sites in that region, progressively added randomly selected human promoters to the list and estimated z. as expected, the z estimate decreased with the addition of random sequences until it stabilized after the addition of a few thousand promoters .
the estimated z's for the upstream  <dig> bp and downstream  <dig> bp regions of all genes in the human genome are  <dig>  and  <dig> , respectively.
effect of the nature and length of training sequences, initial transition probability to the motif  and training of motif profile on trained z
we examined the effect of the nature and the length of training sequences, the initial z  and the effect of the motif profile on trained z. motif profile was kept fixed in this investigation to isolate the effect on z. when trained on tss-n promoters and the initial motif profile, z is inversely proportional to the training promoters' length in the range between 500– <dig> bp.
hence, the quantity z is relatively constant at ~ <dig>  .
it drops slightly between  <dig> to  <dig> bp and then substantially after  <dig> bp due to the lack of κb sites in the shorter training promoters.
when trained on surround-n promoters and the initial motif profile, the trained z is again inversely correlated to the training promoters' length, but the above quantity is higher at ~ <dig>  .
this quantity is even higher at ~ <dig>  when a slightly weaker motif profile corresponding to a surround-n hmm is used  instead of the initial motif profile .
the initial z of tss-n or surround-n hmms in the range between  <dig>  and  <dig>  does not affect the trained z, probably because a global optimum is reached after a few expectation maximization  iterations during training.
varying z with respect to the distance from the tss while scoring
as mentioned in the background, site density decreases in locations that are increasingly further upstream of the tss  <cit> .
also, site density at a location in a sequence is proportional to the transition probability to the motif  at that location.
therefore, while predicting sites in upstream  <dig> bp regions, we modeled z using an exponential functional form such that a sequence close to the tss had a higher site density than a sequence further upstream.
the maximum likelihood estimate of the mean distance of κb sites upstream of the tss was  <dig> bases.
interestingly, the estimate of the mean distance using the median was about the same .
by equating the expressions for the site density per promoter when z had a uniform or exponential function form, the scale factor in the exponential form was z <dig> =  <dig>  for the uniform z =  <dig>  obtained in training .
a different transition probability matrix was generated for each upstream position based on the value of z at that position and was used to calculate occupancy probability.
on the other hand, when identifying κb sites in the downstream  <dig> bp regions, we used location-independent transition probabilities based on the uniform z of  <dig>  obtained in training.
this choice was made due to the paucity of evidence of positional dependence of site density in these regions.
ohmm performs better than a wm
the roc analysis shows that ohmm performs better than a weight matrix  .
while both the hmm and the wm are highly accurate when identifying strong sites, the hmm is more accurate in identifying weak sites.
the segregation of weak sites from site-like sequences is quite difficult due to degeneracy and provides a crucial test.
in this respect, our model far outperforms the wm.
we believe that this superior performance of ohmm is the result of training the threshold in a principled manner to minimize false positives and false negatives.
the positive examples consist of the  <dig> known human κb sites present in upstream  <dig> bp regions  and the negative examples consist of all 10-mers in the upstream  <dig> bp regions in  <dig> randomly selected human genes as described in the methods section.
leave-one-out cross-validation was performed in this analysis.
many κb sites predicted by ohmm are evolutionarily conserved and regulated after nf-κb over-expression
we predicted κb sites in all gene promoters in the human genome and calculated their occupancy probabilities.
two types of data suggest that they may be functional sites.
first, evolutionary conservation scores of κb sites predicted by ohmm are higher than those of  <dig> 10-tuples randomly selected from human promoters, and κb sites with higher hmm occupancy probability have higher evolutionary conservation scores .
secondly, the chicken genes regulated by over-expressed nf-κb proteins in a microarray experiment  <cit>  and their human orthologs are enriched with κb sites predicted by ohmm .
notably, genes regulated in a higher number of comparisons are more enriched with hmm-predicted sites.
also, ohmm predicted more κb sites per regulated gene among genes predicted to contain κb sites probably because true nf-κb targets contain multiple sites.
interestingly, human orthologs of regulated chicken genes are more enriched with predicted nf-κb targets than the chicken genes themselves probably due to the availability of higher quality sequences for humans.
in this experiment, seven nf-κb proteins from different species were over-expressed in chicken dt <dig> pre-b cell lines, and regulated genes were identified by comparing the expression level for each experimental condition against the control.
quantitative measure of biological significance
with the aim of determining the biological significance of the sets of κb target genes predicted by ohmm at the various thresholds, we used pathway analysis to discover the pathways enriched in these κb target genes .
the sum of the negative logarithm of the p-values of the top  <dig> enriched pathways was used as the measure of biological significance.
it is useful to note that only about 50–70% of the genes in each gene set are available for pathway analysis because the rest of the genes are not adequately annotated.
the numbers in figure  <dig>  however, correspond to the number of genes in the entire gene sets.
the gene sets predicted by the hmm at various occupancy probability thresholds are much more biologically significant than randomly selected genes.
moreover, the biological significance reaches a peak at the occupancy probability threshold of  <dig>  .
this implies that the gene sets corresponding to the thresholds greater than  <dig>  have many false negatives, because of which these gene sets do not have enough key target genes to attain high significance.
on the contrary, the gene sets corresponding to the thresholds less than  <dig>  have many more false positives which dilute these sets and lower their biological significance.
thus, the hmm appears to identify sites most accurately in a short window around the threshold of  <dig> .
this observation ties excellently with our justification for training the hmm threshold and the use of occupancy probability as the discriminant function.
correlation with gel shift experiments
a common way to measure affinity of a tf for a particular dna sequence is the gel shift assay, some times also called electrophoretic mobility shift assay or emsa.
in this assay, oligonucleotides  with the particular dna sequence is mixed with the tf at a certain concentration.
at equilibrium, a certain fraction of these oligonulcotides will be bound to the tf.
the oligonucleotides are then passed through a polymer gel under electric field.
dna bound to tf moves at a different rate  than the unbound ones, producing a shift between the two bands.
the relative amount of dna in each of these two bands could be measured, usually by scanning for radioactivity, allowing us to calculate the occupancy probability.
we performed gel shift experiments using double-stranded radio-labeled oligonucleotide probes containing 10-mers derived from several chicken promoters to determine if ohmm predicted the occupancy probability accurately .
we performed the experiments for two nf-κb family members: rela and c-rel.
the transition probability to the motif  in a gel shift experiment is higher than in the cellular context due to higher protein concentration in the gel .
however, because it is difficult to know the protein concentration in the gel, we estimated z for both rela and c-rel .
the estimated z is  <dig> .
occupancy probabilities predicted by the hmm correlate well with the binding affinities of rela and c-rel proteins to these sites at this z .
effect of z on hmm-predicted occupancy probability
we plotted the hmm-predicted occupancy probability with respect to z while keeping the same motif profile .
three characteristics of the dependence between z and occupancy probability stand out:  occupancy probability increases sigmoidally and then saturates as z increases.
occupancy probability of a stronger site  saturates at lower z, and therefore occupancy probability of the stronger site is greater at a particular z.  occupancy probability is influenced by surrounding sequences due to the formation of spurious sites .
hmm reveals the importance of selection of padding sequences of self-overlapping motifs
the occupancy probabilities calculated using the sites in their native chicken promoters did not correlate as well with the experimental binding affinities as those calculated in the previous section using the sites and their padding sequences in the gel shift construct.
when we observed that the difference was due to a c in the padding sequence 3' of the predicted κb sites in the oligonucleotides used for gel shift, we performed a systematic combinatorial analysis using hmm to determine the padding sequences that did not form spurious overlapping binding sites and hence affected native binding the least.
we found that the use of the padding sequences in the above experiment  may misrepresent nf-κb binding.
the predicted occupancy probabilities suggest that a gel shift oligonucleotide containing an a 5' to the 10-mer and a t 3' to the 10-mer in the padding sequence has the least chance of forming spurious binding sites .
any other nucleotide forms a spurious site shifted one position from the 10-mer, and in some cases may even change the binding occupancy by more than 50%.
any of the c, g or t in the 5' end creates a site beginning with cgg..., ggg... or tgg.... an a, c or g in the 3' end creates a site on the opposite strand beginning with tgg..., ggg... or cgg.... in addition, a c at the 3'end also creates a site on the same strand ending with ...ccc.
identification of nf-κb targets in the human genome
as mentioned above, we predicted κb sites in the upstream  <dig> bp and downstream  <dig> bp regions  of all genes in the human genome using the exponential and uniform z respectively .
genes containing at least one κb site with predicted occupancy probability greater than or equal to  <dig>  are designated as putative direct targets of nf-κb.
we also identified cellular pathways, biological functions and diseases in which our predicted nf-κb targets were over-represented.
as expected, we found many known nf-κb target genes with roles in b or t cell receptor signaling, nf-κb signaling, cytokine and chemokine signaling, antigen presentation, acute phase response, or in death receptor and apoptosis signaling among others .
importantly, ohmm also pinpointed several novel candidate nf-κb targets in these and other pathways that have not yet been described to be regulated by nf-κb  <cit> .
this is exemplified by identification of κb sites for diablo, which potentiates some forms of apoptosis, and for the traf family-associated nf-κb activator tank that inhibits nf-κb activation and enhances apoptosis by activating cell death signals and inhibiting survival signals.
although it remains to be determined if these and other genes identified using ohmm are genuine nf-κb transcriptional targets, if confirmed, these could yield important new insights into the less well-characterized but nonetheless important pro-apoptotic activity of nf-κb that has been observed in response to certain stimuli and in certain cells  <cit> .
selected cellular pathways, biological functions and diseases in which our predicted nf-κb targets were over-represented are shown.
the associated predicted nf-κb targets are represented by official human gene symbols.
genes containing κb sites with predicted occupancy probability greater than  <dig>  were used in this analysis.
please see additional file  <dig> for the complete list.
genes known in the literature to be regulated by nf-κb   <cit>  are denoted with *.
similarly, ohmm pinpointed κb sites for ikbkb , which is the key nf-κb activating kinase in the canonical nf-κb signaling pathway , as well as for toll-like receptor  <dig>  that serves as a coreceptor for rna-associated autoantigens and participates in the innate immune response to microbial agents, potentially extending the list of nf-κb-regulated tlrs .
if validated as genuine nf-κb targets, identification of ikkbeta and tlr <dig> by ohmm may have uncovered novel positive feedback loops for amplification of nf-κb signaling.
our analysis suggests that nf-κb could play new roles in some of the pathways in which it is already known to participate, as seen in the protein ubiquitination pathway that can either regulate protein activation/function or target proteins for proteolytic degradation via the proteasome.
for instance, while it is known that nf-κb activates the expression of deubiquitinating enzymes cyld and a <dig>  that negatively regulate nf-κb signaling  <cit>  as well as expression of the proteasome subunit lmp <dig>  and proteasome activators pa28-alpha and pa28-beta  <cit> , ohmm identified the e <dig> ubiquitin-conjugating enzymes ube2h, ube2d <dig> and ube2m that promote protein ubiquitination, as putative nf-κb targets.
incidentally, ube2m was previously observed in microarray studies to be induced by epstein-barr virus that leads to increased nf-κb activity, although it remains to be determined if ube2m is directly controlled by nf-κb  <cit> .
interestingly, our analysis also pinpointed κb sites in btrc .
btrc is the substrate recognition component of a skp1-cul1-f- box protein  e <dig> ubiquitin ligase complex that mediates ubiquitination and subsequent proteasomal degradation of target proteins including the nf-κb inhibitor proteins ikappab alpha , ikappab beta  and ikappab epsilon , allowing free nf-κb dimers to translocate to the nucleus and activate transcription.
btrc also mediates ubiquitination and subsequent proteasomal processing of the phosphorylated nfkb1/p <dig> and nfkb2/p <dig> precursor proteins that respectively participate in the canonical and non-canonical nf-κb signaling pathways  <cit> .
if btrc is a true target of nf-κb, this would suggest that upregulation of btrc by nf-κb could set up a positive feedback loop to amplify degradation of ikappab alpha, thereby providing a novel means for magnifying nf-κb signaling.
overall, these results highlight a potential new role for nf-κb in promoting protein ubiquitination by modulating expression of e <dig> ubiquitin-conjugating enzymes and e <dig> ubiquitin ligases, thus expanding on its previously documented role in inducing expression of deubiquitinases, proteasome activators and proteasome subunits.
ohmm also revealed new insights into the cross-talk between nf-κb and other signaling pathways.
for example, we identified κb sites in several key mediators in the notch signaling pathway that is involved in cell-cell communications to regulate a broad spectrum of cell-fate determinations.
these include the delta-like  <dig> ligand for notch receptors dll <dig>  the notch <dig> receptor, transcriptional regulator rbp that acts as a transcriptional repressor in absence of notch but is a transcriptional activator when associated with activated notch, and mastermind-like  <dig>  that serves as a transcriptional coactivator for notch.
prior studies have suggested multiple cross-talk mechanisms between the notch and nf-κb signaling pathways, most of them suggesting activation of nf-κb downstream of notch, as well as protein-protein interactions between components of these two pathways .
our finding that key mediators of notch signaling harbor κb sites supports the notion that nf-κb might in turn modulate notch signaling to influence cell fate determination during development, immunity and cancer.
this would be consistent with recent work showing that nf-κb stimulates the expression of notch targets hes- <dig> and deltex- <dig> to stimulate marginal zone  b-cell development  <cit>  and with prior observations from our group that nf-κb can trigger notch signaling by inducing expression of the notch ligand jagged <dig>  although it remains to be determined if its effect on jagged <dig> expression is direct  <cit> .
other genes uncovered by ohmm could be intimately associated with nf-κb's role in transcriptional regulation, and in particular, with its recently uncovered role in gene-specific transcriptional repression  <cit> .
for example ohmm predicts κb sites for deacetylases hdac <dig> and sirt <dig>  as well as for transcriptional corepressor sin3a.
while nf-κb is known to engage in protein-protein interactions with histone deacetylases and corepressors to modulate gene transcription, ohmm raises the possibility of a new mode of action for nf-κb in this context.
also highlighted were candidate targets in xenobiotic metabolism, including κb sites in sulfotransferase sult1c <dig>  aldehyde dehydrogenase  <dig> family gene aldh3b <dig>  and transcription factor nuclear factor erythroid derived 2-like  <dig>  that regulates the oxidative stress response, in addition to the previously documented nf-κb target glutathione s-transferase .
this points to a possible role for nf-κb in drug metabolism, multidrug resistance and detoxification of poisonous compounds, with possible impact for treatment of infectious diseases, anti-cancer therapy and/or environmental science.
finally, ohmm pinpointed possible involvement of nf-κb in new pathways, such as the neurotrophin/trk signaling cascade, in which its role is emerging.
in addition to the recently recognized nf-κb target phosphatidylinositol 3-kinase catalytic subunit pik3c <dig>  <cit> , ohmm identified κb sites for its regulatory subunit pik3r <dig>  as well as for the ras-related g-protein rras, ribosomal protein s <dig> kinase  <dig> rps6ka <dig> and transcription factor cyclic amp-responsive element-binding protein  <dig> .
in this regard, creb <dig>  was recently described to be controlled by nf-κb  <cit>  as is the critical neurotrophin nerve growth factor ngf that controls neuronal cell fate  <cit> .
these results may warrant further investigation in view of the emerging role of nf-κb in the regulation of neuronal survival and function in the nervous system, its effect on cognition and behavior and its suggested roles in epilepsy, stroke, alzheimer's, parkinson's diseases and other neurological disorders .
verification of some of the predicted nf-κb sites in the human genome
gel shift experiments were performed on several of the predicted nf-kb sites using two of the nf-κb family members: rela and c-rel.
three low occupancy sequences were chosen as negative control and the known nfkb1a site was taken as the positive control.
we chose the predicted target sites to be tested from some of the top categories that were extracted by pathway analysis.
the sites to be tested were from the promoters of following genes: ikbkb and tlr <dig>  involved in the nf-kb signaling pathway; diablo, tank, bcl2l <dig> and bnip <dig>  involved in the death receptor signaling and apoptosis categories; hdac <dig>  sirt <dig> and sin3a, involved in cell cycle g1/s and transcriptional regulation; and btrc, ube2d <dig> and ube2m, involved in ubiquitination pathway.
all of these sites show stronger gel shift compared to the negative controls.
to understand which whether predicted occupancy is a good indicator of binding compared to, say, evolutionary sequence conservation, we decided to classify these sites as the following as shown in figure  <dig>  the conservation score greater than 90% is considered high, between 90% and 80% is considered to be medium and less 80% is considered to be low.
similarly, sites with occupancy greater than  <dig>  is considered to be high occupancy and the rest are considered to be low occupancy.
the tested sites, organized according to this classification, are listed below.
high conservation and high occupancy: ifbkb, ube2d <dig> and sin3a
high conservation and low occupancy: tank- <dig>  tank- <dig>  bcl2l <dig> and ube2m
medium conservation and high occupancy: sirt1
medium conservation and low occupancy: tlr <dig>  hdac <dig> and btrc
low conservation and high occupancy: bivip <dig> and diablo
in general, the high occupancy sites are better binders, the exceptions being sin3a , hdac <dig>  and tlr <dig> .
for hdac <dig>  the occupancy is only slightly lower than the cut-off threshold.
in comparison, most highly conserved low occupancy sites have weak gel shifts, suggesting ohmm occupancy measure is a better predictor of binding affinity than the sequence conservation score.
discussion
while the use of an hmm is thus advantageous in identifying sites, we have discovered that an hmm offers tremendous benefits for the special case of identifying self-overlapping motifs.
when putative sites overlap, the overall occupancy probability in a position in the dna sequence depends upon the strength of all sites containing that position.
calculation of the overall occupancy probability over the entire sequence using the gamma variables in the hmm method enables us to compute the occupation probability in presence of overlapping sites.
in this regard, our approach has a distinct advantage over the two hmm scoring methods in vogue at the present time, viz.
the likelihood method and the viterbi method.
the likelihood of a sequence scored over a long stretch of sequence is greater when more sites are present in that sequence.
it is perfect for discovering clusters of sites, as evident from the work on fly regulatory modules.
the likelihood score by itself does not indicate the location of sites in a promoter.
in contrast, the viterbi method proclaims the presence of sites at positions where the state path with the highest probability contains the motif state.
it, however, fails to consider overlapping sites because only the best site would be present in the optimal state path.
training ohmm is a non-trivial two-step process.
in the first step, we trained both the emission and transition probabilities using short sequences rich in known sites with the aim of accurately estimating the motif profile.
in the second step, we kept the emission probabilities constant, and trained the transition probabilities on promoters containing known sites as well as random promoters to accurately estimate the transition probability to the motif z reflecting the site density in the promoters of all genes in the human genome.
unlike most hmms in the literature, we successfully trained the hmm emission probabilities without requiring pre-alignment of training sites.
one interesting feature is that trained z is inversely related to the strength of the motif profile.
in other words, when the motif profile is kept constant and only z is trained, the weaker the motif profile used , the higher the trained z. this is probably the result of the compensating mechanism between the motif profile and z discussed above.
interestingly, this competition between z and the motif profile also determines the occupancy probability as we saw before.
the goal of the second training step was to estimate the transition probability to the motif z reflecting the site density in the promoters of all genes in the human genome.
this z corresponds to the appropriate threshold when identifying sites in all human promoters.
obviously, z trained in the first step was not appropriate due to the high site density in the site-rich training sequences.
instead of training z on all human promoters, which is computationally expensive, we train it on a sufficient number of promoters to get a reasonable estimate of z. our procedure of beginning with the human promoters containing known sites as the training set and progressively adding promoters of randomly selected human genes to the training set leads to decreasing z, until the training set reached a few thousand promoters and z stabilized .
we expect that, for a particular transcription factor, the speed of convergence of this procedure depends upon overall number of target genes of the factor: larger the number of targets/sites, quicker the convergence.
one shortcoming of most site identification methods is that they assume equal probability of the site in a window of arbitrary length around the tss.
the majority of known κb sites are located just upstream of the tss in gene promoters and the number of known κb sites decreases further upstream.
specifically, of the  <dig> known κb sites upstream of the tss,  <dig> are located within  <dig> bp and  <dig> are located within  <dig> bp of the tss.
liu et al.
have also made a similar observation in the promoters of nf-κb-regulated immune genes  <cit> .
to counter the claim that such an observation for binding sites may be due experimental bias, we cite tabach et al.
who showed in a wide-scale bioinformatic study that functional binding sites are more likely to be present in the  <dig> bp region upstream of the tss than any other upstream region for most human transcription factors and specifically for nf-κb  <cit> .
they defined functional binding sites as those over-represented in functionally related genes  and conserved in related species.
to bolster their conclusion, they showed location dependence of binding sites for the transcription factor myocardin in a controlled experiment.
xie et al.
also arrived at a similar conclusion based on binding site conservation  <cit> .
even though the exact reasons for the occurrence of such a phenomenon are not known at the present time, a better interaction of the transcription factor with the transcription machinery if it is bound close to the tss and the low density of nucleosomes near the tss  <cit>  heuristically explain why this phenomenon may occur.
while we agree that a great number experiments need to be conducted to definitively prove the location dependence of functional κb sites, we feel that a site identification method must be able to take it into account.
noticing that site density decreases sharply with the upstream disease from the tss and that transition probability to the motif z is proportional to site density, we modelled z using an exponential functional form such that the probability of identifying a site decreases further upstream but never reaches zero .
even though this approach may fail to identify sites in distal promoters and enhancers, we believe that it allows site search in large upstream regions without identifying too many false positives.
a κb site exerts influence on the occupancy probabilities at the positions surrounding the site in either direction.
because overlapping binding sites are usually present when the scoring window is moved by one position, occupancy probability stays high at all positions in that window.
when the window is moved by more positions, the occupancy probability at the new positions dips slightly below the average due to the high motif probability at the site.
the occupancy probability returns to the average background value when the scoring window moves by ten positions.
the relationship between the occupancy probabilities of sites located close to each other is quite instructive.
when two sites are in tandem without any space between them, occupancy probabilities of both of them are lower due to the small motif-to-motif transition probability.
occupancy probabilities are not very high even when the sites are one position apart because a window shift by one position from a κb site usually contains an overlapping κb site.
the sites need to be at least two positions apart so that they do not exert significant influence on each other's occupancy probabilities.
please note that in any case, the overall occupancy in a region containing two nearby sites is quite high.
a strong correlation exists between the occupancy probabilities predicted by ohmm and the in vitro binding affinity of the nf-κb family members c-rel and rela for oligonucleotides in a gel shift experiment .
this validates our physical binding model because the hmm-predicted occupancy probabilities appear to correspond to observed binding affinities.
in the process of the study, we found how the padding sequence affects the results.
the current practice for performing gel shift experiments in the nf-κb community consists of using particular padding sequences around the 10-mer corresponding to a potential κb site .
the commonly used padding sequences, however, can form spurious κb sites due to the self-overlapping nature of κb motif, and hence the experiment may not represent binding of the 10-mer in the native promoter.
ohmm offers guidance on the selection of correct padding sequences when designing experiments.
a gel shift oligonucleotide should have an a 5' to the κb site and a t 3' to the κb site for minimum interference, and that a 3' c should be avoided at all cost.
ideally, oligonucleotides containing a few bases corresponding to those surrounding the κb site in the promoter of the gene should be used to pad the site.
this will capture the effects of all the overlapping putative binding sites in the native promoter and will avoid creation of artefacts based on the nucleotides present in the padding sequence.
the nf-κb family of transcription factors consists of some of the most important proteins for an organism's development and survival.
in resting cells, inactive nf-κb dimers reside in the cytoplasm, bound to inhibitors in the iκb family  <cit> .
many different stimuli can activate latent nf-κb/iκb complexes, including pro-inflammatory cytokines that trigger rapid activation of the classical  nf-κb signaling cascade via the ikk kinase complex  <cit> .
consequent phosphorylation and degradation of iκb via the proteasome culminates in nuclear translocation of active nf-κb dimers and their binding to κb sites.
this commonly results in transcriptional activation, or occasionally in repression, of genes important for the immune and inflammatory response, cell proliferation, adhesion, angiogenesis and inhibition of apoptosis  <cit> .
the non-canonical  nf-κb signaling cascade is activated via the lymphotoxin and baff receptors, leading to activation of ikkα, nuclear translocation of p52/relb complexes and plays a role in lymphoid organogenesis and b-cell function.
hundreds of nf-κb-regulated genes have been identified .
among them, binding of activated nf-κb to the iκbα promoter can trigger its resynthesis, giving rise to an autoregulatory loop that terminates the activation process  <cit> .
consequently, under normal conditions, rel/nf-κb activation is tightly regulated and transient.
however sustained activation of rel/nf-κb can promote cell survival, proliferation and/or cell death depending upon the nature and the extent of the stimulus, and is implicated in many pathological conditions including immune system disorders, chronic inflammation and cancer  <cit> .
although several nf-κb target genes have been identified, whether other nf-κb-regulated genes are its direct or indirect targets is not known.
in addition to identifying many of these known nf-κb target genes, ohmm pinpointed several novel candidate nf-κb targets that are yet to be shown as controlled by nf-κb.
while further studies will be needed to determine which of the candidate genes identified by ohmm are genuine transcriptional targets of nf-κb, these are likely to yield important new insights into the less well-characterized roles of nf-κb in the induction of apoptotic cell death, its recently uncovered role in transcriptional repression, the regulation of protein ubiquitination and its role in other signaling pathways, to name a few.
in addition to identifying multiple genes with key roles in haematological and immune cell development and function, ohmm pinpointed κb sites in genes with roles in drug metabolism, skeletal and muscle system development and function, as well as in disease conditions ranging from immune and inflammatory disorders to infectious diseases, cancer, skeletal muscular disorders and neurological diseases .
these findings are consistent with experimental evidence suggesting a role for nf-κb in immunological and inflammatory diseases, cancer and therapy-resistance, in skeletal myogenesis and cachexia, as well as in cognition, behavior and neurological disorders like epilepsy, stroke, alzheimer's, parkinson's diseases, and amyotrophic lateral sclerosis  .
ohmm is thus a unique tool to shed new light on the multifaceted biological functions of nf-κb and the mechanisms involved, and to identify critical mediators of dysregulated nf-κb activity in human disorders that may ultimately help to develop new approaches for treatment.
let us finish the discussion by mentioning possible shortcomings of the method.
ohmm requires a complicated two-step training procedure, which might only work practically for factors with many targets.
we have not considered insertions or deletions inside the motif, making ohmm well suited for the simultaneous training of emission and transition probabilities.
in addition, we have considered the transition probability to the motif z in the upstream region to decrease exponentially during the scoring procedure.
this density function may not be accurate, especially for enhancers.
we will need many more sites to establish a more accurate density function.
CONCLUSIONS
our method is the first attempt to deal specifically with self-overlapping binding motifs, providing guidance in the selection of padding sequences in gel shift experiments.
when the hmm is considered as a biophysical model, it interprets the weight matrix as binding energy, the transition probability to the motif as transcription factor concentration and the gamma variable as occupancy probability.
because it uses occupancy probability as the discriminant function, it learns the associated natural threshold in a principled manner during the training procedure.
another unique feature of ohmm is that transition probabilities change with sequence location to reflect site density.
transcription factor's occupancy probability on a site can be accurately estimated using our method, as seen by the high correlation with experimental binding affinities.
high evolutionary conservation scores and enrichment in experimentally regulated genes suggest that κb sites predicted by our method might be functional.
our results may provide important new insights into the function and regulation of nf-κb and uncover possible new biological roles for this important transcription factor family.

 identification of short, degenerate dna sequences that bind to a transcription factor is a difficult problem. the particular sub-problems of identification of sites corresponding to self-overlapping motifs, determination of a threshold for this purpose, biophysical interpretation of statistical quantities used in probabilistic site identification methods and estimation of occupancy of these sites by a transcription factor have not been addressed in detail before. these problems form the focus of this study. several transcription factors bind to self-overlapping dna motifs. when a motif overlaps with itself with a particular shift, the transcription factor can bind in more than one sequence windows, i.e. to overlapping sites. examples include drosophila developmental transcription factor hunchback, worm pha-, human sp- c/ebpalpha, yeast adr mig chicken cdx- arabidopsis agamous, etc.. furthermore, when binding by the transcription factor in either orientation is permissible, the corresponding dna site and its reverse complement can be considered as two different types of sites that are overlapping, as long as the motif is not exactly palindromic. problems similar to the identification of overlapping sites also arise in the context of prediction of nucleosome positioning. self-overlapping motif of nf-κb family of factors the nf-κb family of transcription factors is a prominent example of transcription factors with self-overlapping motifs. these factors are key mediators of the cellular response to infection, injury, inflammation, or other stress conditions that lead to rapid alterations in cellular gene expression. the vertebrate members include the c-rel, rela, relb, nfkb1/p and nfkb2/p proteins that form homo- or heterodimers with one another and bind to similar sites. these nf-κb binding sites, with the consensus being gggrnnyycc, often overlap because a κb site contains multiple g's at the 5' end and multiple c's at the 3' end with a high probability. for a good κb site, when the sequence window is shifted by one base in either the 5' or the 3' direction the resulting sequence is also a putative κb site. moreover, we can have additional possibility of self-overlap because the reverse complement of a κb site is often a κb site, allowing functional binding in the opposite direction. many computational methods do not take into account these alternative binding modes while scoring a candidate regulatory site. ignoring the self-overlapping nature of the motif can also complicate the interpretation of in vitro experiments. for example, a 3' padding sequence, starting with nucleotide c, in a gel shift experiment can form a spurious strong κb site and thus confer higher binding affinity to the experimental sequence even when the test sequence, in the context of the native promoter, makes for a very weak κb site. a computational method that naturally considers multiple binding windows and help in the design of padding sequences is useful for avoiding this problem. the commonly used site identification methods often assume equal probability of a site everywhere in a particular window around the transcription start site and thus fail to adjust according to the observed distribution of sites within the gene structure. whereas proximal promoters, up to bp upstream of the tss, contain several sites, further upstream distal promoters contain fewer sites. while identifying sites, a method should adjust its parameters to reflect this change in the density of sites with respect to the distance from the tss. a computational method that links to the biophysical process of transcriptional regulation can offer more biological insight than most of the current site identification methods. such a method can interpret statistical quantities as biophysical variables like binding energy and transcription factor concentration. it can also estimate how often a transcription factor is bound to the gene promoter, i.e. its occupancy probability. interestingly, the use of occupancy probability as a discriminant function during site identification has a number of advantages. one advantage of this approach is that it goes beyond just site identification and provides a measure of the importance of a site in modulating gene regulation. for example, for an activator, highly occupied sites may result in high gene expression and thus may be biologically more significant. most methods that identify transcription factor binding sites or regulatory modules, however, leave unexplored the occupancy of the transcription factor on the promoter and its impact on gene regulation. the lack of biophysical interpretation prevents these methods from using a natural threshold for classifying sequences into sites. they declare that a sequence is a site if its score is above an arbitrary threshold or is statistically significant compared to the score of the random background. they often have a high threshold resulting in several false negatives or a low threshold resulting in several false positives. thus, identification of weak sites is difficult. on the contrary, saturation effects in the physical occupancy provides a natural threshold and this fact has been utilized to develop class support vector machines like qpmeme. alternatively, hidden markov models, when interpreted as biophysical models, overcome this problem by learning the threshold based on the transition probability to the motif state, namely the hidden state associated with a sequence pattern which is distinct from the background. hmm for binding site identification hmms, now popular in bioinformatics for more than a decade, have been used in two different ways for motif identification. for identification of one or more occurrences of non-overlapping sites:'profile hmms', originally designed to model protein domains, have more recently been used to identify binding sites of transcription factors, for example, of camp receptor protein in cyanobacterium anabaena, liver x receptor and creb. a profile hmm library was built using transfac sequences to classify transcription factors. in a profile hmm, each position within a motif has three states. a match state is associated with a nucleotide being present at that position and has corresponding emission probabilities. a deletion state corresponds to absence of any nucleotide at that position. an insertion state allows for insertion of nucleotides between the current position and the next position within the motif, and has its own emission probabilities. for identification of cis-regulatory modules that contain multiple sites of different types: this is usually performed using'motif hmms'. in a motif hmm, the entire motif, i.e. all positions within the motif, is represented by one state. different states correspond to different motif types. phylogenetic conservation has been incorporated in motif hmms to reduce false positives. both types of hmms, however, do not generally train their parameters properly. a profile hmm has a complicated architecture and requires a large number of parameters as a consequence. because the number of known sites is generally small, training of a profile hmm using known sites in their native promoters is usually not possible. therefore, a profile hmm is generally trained using pre-aligned sites. because the whole promoters containing the training sites are not used, transition probability to the motif z from background is not trained, and the relationship between transition and emission probabilities is likewise not captured. motif hmms, on the other hand, are more focused on identifying motifs of multiple types. but while they attempt to estimate the transition probability to the motif z, they generally use only the promoters containing the motifs for training z and thus overestimate z. moreover, they usually train emission probabilities of motifs separately using pre-aligned training sites, thus ignoring the effect of z on emission probabilities. many of these hmms use the likelihood or the viterbi algorithm for scoring, and thus end up using an arbitrary classification threshold. furthermore, they leave their relationship with biophysical models rather obscure and thus fail to calculate the occupancy probability. these shortcomings are in addition to their basic failure to explicitly consider the effect of multiple overlapping sites in the viterbi method. methods specifically developed for κb site identification in an original approach to specifically identify κb sites, udalova et al. developed a principal coordinate model. they determined relative binding affinities of κb sites using experimental quantitative binding data. they selected a subset of the possible variants of the fully palindromic nf-κb binding consensus sequence ggrrnnyycc such that no variant differed from the selected sequences or their reverse complements by more than one nucleotide. they mapped these sequences to a euclidean space and used the largest principal components for least-square linear regression of the logarithm of binding affinity in a gel shift assay or microarray. this model automatically incorporated effects of interactions between base pair positions in the binding site, its predictions were highly correlated with experimental binding data and it identified site positions responsible for differential binding of nf-κb family members consistent with crystallographic studies. however, the model's disadvantages are that it requires experimental quantitative binding data of all selected sequences and only includes variants of the consensus sequence. because several known κb sites do not fit the consensus sequence, inclusion of all possible 10-mer variants for this model will require binding experiments with a large number of sequences, making this model infeasible. furthermore, this model also suffers from the limitations of site identifications described above. alternative interpretation of hmm as a physical binding model even though in the context of site identification an hmm is usually interpreted as a sequence generative model, we focus on a less familiar interpretation as a physical binding model of a transcription factor on dna. this interpretation leads us to transform the statistical hmm model into a biophysical one. we can then determine the occupancy probability of a transcription factor on a dna sequence, and think of the prior probability of the motif as a measure of transcription factor concentration and the weight matrix as a measure of binding energies. more importantly, the biophysical model offers a principled threshold for classifying sequences into sites. the statistical mechanics model was utilized in while the equivalent hmm representation is used in. although we do not know of a paper explicitly describing the equivalence of these two approaches, the work on nucleosome positioning by segal et al. utilizes this relationship. an hmm is commonly used as a generative model of a spatial or time series sequence in the machine learning field. in our context, the hmm consisting of two kinds of states – the background state and motif states – generates a sequence from 5' to 3'. at any position in the sequence, the hmm determines the probabilities of the motif and background states at the previous position, calculates the probability of either state at the current position using the transition probabilities and generates the new nucleotide based upon the hidden states' emission probabilities. we could also consider the hmm as a physical binding model that estimates the occupancy probability of a transcription factor on a particular position of a dna sequence, i.e. how often the transcription factor is bound to that position of the dna sequence. the relation between the two interpretations of an hmm, are as follows. in our particular hmm, the start positions of a motif state, with a particular orientation, specifies the hidden states. one could ask define the conditional probability distribution of the hidden states, or, equivalently these positions, is given the parameters and the observed sequence. it turns out that this conditional distribution is exactly the same as that of the edge positions of bound transcription factors on dna in the statistical mechanical model. the occupancy probability at a position is the probability of finding a motif state at that position in the hmm language. if we ignore the possibility of binding in multiple overlapping windows, then, the probability at the jth position of sequence s being the beginning of a motif is, where wj is the weight matrix score of the motif starting at the jth position of the sequence. we claim that this quantity can be mapped to the occupancy probability, as we will see below. thus, the two factors determining occupancy probability are the transition probability to the motif z and the measure of distinctness of the emission probabilities of the motif from that of the background. both these factors need to be high for the transcription factor to be bound to a particular position in a dna sequence with high probability. for example, even if the weight matrix score is high, occupancy probability cannot reach one if z is really small. a site identification method based only on a weight matrix has no way of dealing with this interplay with z. the hmm training techniques offer two advantages over the calculations made using first principles. first, the hmm baum-welch procedure trains the transition and emission probabilities. this produces optimized values of z and the weight matrix, which are essential for an accurate estimation of occupancy probability. secondly, as mentioned above, the hmm training procedure also calculates an intermediate variable called the gamma variable. the occupancy probability at a particular position is given simply by the gamma variable of the motif state at that position:, where and are the gamma values of the motif and background states at that position, respectively. thus, the calculation of occupancy probability in a new sequence requires a simple extension of the scoring procedure in which the gamma variable is computed just like during the training procedure. in a biophysical model, based on thermodynamics principles, occupancy probability can be written as, where is the concentration of a free transcription factor at equilibrium, ej is the binding energy at position j, where kb is the boltzmann's constant, and t is the absolute temperature. a comparison of this equation with the equation of probability obtained from the hmm gives us wj = - βej and z =. thus, the weight matrix represents binding energy. in addition, the transition probability to the motif z corresponds to the free transcription factor concentration. as the transcription factor concentration increases, the transition probability to the motif increases and we expect higher occupancy by the transcription factor on the dna. one main focus of this study is to determine a good threshold for the classification of sequences into binding sites. use of an arbitrary threshold while using a purely statistical quantity as a discriminant function is prone to a high number of false positives or false negatives. a good discriminant function has a physical interpretation, it relates to the biological significance of a sequence and it offers a natural threshold for classification. moreover, the classification method needs to be able define the above threshold based on training sequences. we use occupancy probability as a discriminant function to classify sequences into sites because it has the following benefits. first, unlike a purely statistical entity, it has a straightforward physical interpretation. a transcription factor's occupancy on the promoter determines gene expression. second, highly occupied binding sites may be physiologically more significant. thus, occupancy probability not only helps in classifying sequences as binding sites but also offers insight into their influence on gene expression. third, occupancy probability has a "natural" threshold at due to its fermi-dirac distribution with respect to binding energy. sequences with binding energy less than the chemical potential have occupancy probability greater than and hence can be classified as sites. because of our focus on occupancy, we call this model ohmm. the greatest benefit of the transformation of an hmm to a biophysical model is that it enables the hmm to learn the threshold for classification of sequences into sites in a principled way. we have seen above that the occupancy probability depends upon the transition probability to the motif z and the weight matrix, both of which are trained by an hmm. thus, when an hmm uses occupancy probability as a discriminant function, it learns the natural threshold based on the training sequences. because of this accurate estimation of the threshold, an hmm is expected to identify weak sites much more accurately with a fewer false positives. in summary, we describe an hmm, applied to nf-κb, that identifies overlapping sites, learns the threshold in a principled manner by training emission probabilities using known sites in their native promoters and training transition probabilities using promoters in an entire genome, estimates the transcription factor's occupancy probability, and adjusts parameters reflecting the distance from the tss while scoring. the structure of hidden states in ohmm is indicated in figure we demonstrate that the occupancy probability predicted by ohmm correlates well with the binding affinity in gel shift experiments. high evolutionary conservation scores and enrichment in experimentally regulated genes suggest that the predicted sites might be functional.  hmm training training ohmm consists of estimating two sets of parameters: the emission probabilities of each motif state and the background, and the transition probabilities that depend upon the transition probability to the motif. while an hmm needs to be trained using site-rich sequences to learn the motif profile, training on random sequences is required to learn z reflecting the site density in the promoters of all genes in the human genome. we therefore estimated hmm parameters in two steps. we first trained all hmm parameters using site-rich sequences and thus learned the emission probabilities. we then fixed the emission probabilities and trained the transition probabilities using randomly selected human promoters to estimate z separately for the tss-800:tss and tss:tss+ regions, where tss is the transcription start site. the reason for separate z estimation in these two types of regions is that κb site density is different upstream and downstream of the tss. training the emission probabilities in the first training step, we used two types of site-rich sequences of different lengths and various initial z's to train all hmm parameters and determine the emission probabilities to be used in further analysis. trained motif profiles of "tss-n hmms," i.e. hmms trained on promoters consisting of n bases upstream of the tss of human genes known to contain a κb site, appear similar to the background emission probabilities regardless of the promoter length and initial z. on the other hand, trained motif profiles of "surround-50" or "surround- hmms," i.e. hmms trained on or bp sequences each consisting of a known κb site and surrounding region, with a reasonable initial z, are distinct from the background. they are also distinct from the initial motif profile, as their symmetrical kullback-leibler divergences with respect to the initial motif profile are high. the trained motif profiles are slightly weaker than the initial motif profile, i.e. more similar to the background. in a surround- or surround- hmm, any initial z between and results in the same trained motif profile, indicating that perhaps a local optimum is reached. trained motif profiles of surround- hmms, however, appear more and more like the background as the initial z increases above. trained motif profiles of surround- hmms appear similar to the background regardless of the initial z. we used the trained motif profile of the surround- hmm for further analysis. estimation of the transition probability to the motif in the second training step, we estimated z by training the transition probabilities of the above surround- hmm using random sequences while keeping its emission probabilities constant. we estimated z separately for the upstream bp and downstream bp regions. in each case, we determined human promoters with known κb sites in that region, progressively added randomly selected human promoters to the list and estimated z. as expected, the z estimate decreased with the addition of random sequences until it stabilized after the addition of a few thousand promoters. the estimated z's for the upstream bp and downstream bp regions of all genes in the human genome are and, respectively. effect of the nature and length of training sequences, initial transition probability to the motif and training of motif profile on trained z we examined the effect of the nature and the length of training sequences, the initial z and the effect of the motif profile on trained z. motif profile was kept fixed in this investigation to isolate the effect on z. when trained on tss-n promoters and the initial motif profile, z is inversely proportional to the training promoters' length in the range between 500– bp. hence, the quantity z is relatively constant at ~. it drops slightly between to bp and then substantially after bp due to the lack of κb sites in the shorter training promoters. when trained on surround-n promoters and the initial motif profile, the trained z is again inversely correlated to the training promoters' length, but the above quantity is higher at ~. this quantity is even higher at ~ when a slightly weaker motif profile corresponding to a surround-n hmm is used instead of the initial motif profile. the initial z of tss-n or surround-n hmms in the range between and does not affect the trained z, probably because a global optimum is reached after a few expectation maximization iterations during training. varying z with respect to the distance from the tss while scoring as mentioned in the background, site density decreases in locations that are increasingly further upstream of the tss. also, site density at a location in a sequence is proportional to the transition probability to the motif at that location. therefore, while predicting sites in upstream bp regions, we modeled z using an exponential functional form such that a sequence close to the tss had a higher site density than a sequence further upstream. the maximum likelihood estimate of the mean distance of κb sites upstream of the tss was bases. interestingly, the estimate of the mean distance using the median was about the same. by equating the expressions for the site density per promoter when z had a uniform or exponential function form, the scale factor in the exponential form was z = for the uniform z = obtained in training. a different transition probability matrix was generated for each upstream position based on the value of z at that position and was used to calculate occupancy probability. on the other hand, when identifying κb sites in the downstream bp regions, we used location-independent transition probabilities based on the uniform z of obtained in training. this choice was made due to the paucity of evidence of positional dependence of site density in these regions. ohmm performs better than a wm the roc analysis shows that ohmm performs better than a weight matrix. while both the hmm and the wm are highly accurate when identifying strong sites, the hmm is more accurate in identifying weak sites. the segregation of weak sites from site-like sequences is quite difficult due to degeneracy and provides a crucial test. in this respect, our model far outperforms the wm. we believe that this superior performance of ohmm is the result of training the threshold in a principled manner to minimize false positives and false negatives. the positive examples consist of the known human κb sites present in upstream bp regions and the negative examples consist of all 10-mers in the upstream bp regions in randomly selected human genes as described in the methods section. leave-one-out cross-validation was performed in this analysis. many κb sites predicted by ohmm are evolutionarily conserved and regulated after nf-κb over-expression we predicted κb sites in all gene promoters in the human genome and calculated their occupancy probabilities. two types of data suggest that they may be functional sites. first, evolutionary conservation scores of κb sites predicted by ohmm are higher than those of 10-tuples randomly selected from human promoters, and κb sites with higher hmm occupancy probability have higher evolutionary conservation scores. secondly, the chicken genes regulated by over-expressed nf-κb proteins in a microarray experiment and their human orthologs are enriched with κb sites predicted by ohmm. notably, genes regulated in a higher number of comparisons are more enriched with hmm-predicted sites. also, ohmm predicted more κb sites per regulated gene among genes predicted to contain κb sites probably because true nf-κb targets contain multiple sites. interestingly, human orthologs of regulated chicken genes are more enriched with predicted nf-κb targets than the chicken genes themselves probably due to the availability of higher quality sequences for humans. in this experiment, seven nf-κb proteins from different species were over-expressed in chicken dt pre-b cell lines, and regulated genes were identified by comparing the expression level for each experimental condition against the control. quantitative measure of biological significance with the aim of determining the biological significance of the sets of κb target genes predicted by ohmm at the various thresholds, we used pathway analysis to discover the pathways enriched in these κb target genes. the sum of the negative logarithm of the p-values of the top enriched pathways was used as the measure of biological significance. it is useful to note that only about 50–70% of the genes in each gene set are available for pathway analysis because the rest of the genes are not adequately annotated. the numbers in figure however, correspond to the number of genes in the entire gene sets. the gene sets predicted by the hmm at various occupancy probability thresholds are much more biologically significant than randomly selected genes. moreover, the biological significance reaches a peak at the occupancy probability threshold of. this implies that the gene sets corresponding to the thresholds greater than have many false negatives, because of which these gene sets do not have enough key target genes to attain high significance. on the contrary, the gene sets corresponding to the thresholds less than have many more false positives which dilute these sets and lower their biological significance. thus, the hmm appears to identify sites most accurately in a short window around the threshold of. this observation ties excellently with our justification for training the hmm threshold and the use of occupancy probability as the discriminant function. correlation with gel shift experiments a common way to measure affinity of a tf for a particular dna sequence is the gel shift assay, some times also called electrophoretic mobility shift assay or emsa. in this assay, oligonucleotides with the particular dna sequence is mixed with the tf at a certain concentration. at equilibrium, a certain fraction of these oligonulcotides will be bound to the tf. the oligonucleotides are then passed through a polymer gel under electric field. dna bound to tf moves at a different rate than the unbound ones, producing a shift between the two bands. the relative amount of dna in each of these two bands could be measured, usually by scanning for radioactivity, allowing us to calculate the occupancy probability. we performed gel shift experiments using double-stranded radio-labeled oligonucleotide probes containing 10-mers derived from several chicken promoters to determine if ohmm predicted the occupancy probability accurately. we performed the experiments for two nf-κb family members: rela and c-rel. the transition probability to the motif in a gel shift experiment is higher than in the cellular context due to higher protein concentration in the gel. however, because it is difficult to know the protein concentration in the gel, we estimated z for both rela and c-rel. the estimated z is. occupancy probabilities predicted by the hmm correlate well with the binding affinities of rela and c-rel proteins to these sites at this z. effect of z on hmm-predicted occupancy probability we plotted the hmm-predicted occupancy probability with respect to z while keeping the same motif profile. three characteristics of the dependence between z and occupancy probability stand out: occupancy probability increases sigmoidally and then saturates as z increases. occupancy probability of a stronger site saturates at lower z, and therefore occupancy probability of the stronger site is greater at a particular z. occupancy probability is influenced by surrounding sequences due to the formation of spurious sites. hmm reveals the importance of selection of padding sequences of self-overlapping motifs the occupancy probabilities calculated using the sites in their native chicken promoters did not correlate as well with the experimental binding affinities as those calculated in the previous section using the sites and their padding sequences in the gel shift construct. when we observed that the difference was due to a c in the padding sequence 3' of the predicted κb sites in the oligonucleotides used for gel shift, we performed a systematic combinatorial analysis using hmm to determine the padding sequences that did not form spurious overlapping binding sites and hence affected native binding the least. we found that the use of the padding sequences in the above experiment may misrepresent nf-κb binding. the predicted occupancy probabilities suggest that a gel shift oligonucleotide containing an a 5' to the 10-mer and a t 3' to the 10-mer in the padding sequence has the least chance of forming spurious binding sites. any other nucleotide forms a spurious site shifted one position from the 10-mer, and in some cases may even change the binding occupancy by more than 50%. any of the c, g or t in the 5' end creates a site beginning with cgg..., ggg... or tgg.... an a, c or g in the 3' end creates a site on the opposite strand beginning with tgg..., ggg... or cgg.... in addition, a c at the 3'end also creates a site on the same strand ending with...ccc. identification of nf-κb targets in the human genome as mentioned above, we predicted κb sites in the upstream bp and downstream bp regions of all genes in the human genome using the exponential and uniform z respectively. genes containing at least one κb site with predicted occupancy probability greater than or equal to are designated as putative direct targets of nf-κb. we also identified cellular pathways, biological functions and diseases in which our predicted nf-κb targets were over-represented. as expected, we found many known nf-κb target genes with roles in b or t cell receptor signaling, nf-κb signaling, cytokine and chemokine signaling, antigen presentation, acute phase response, or in death receptor and apoptosis signaling among others. importantly, ohmm also pinpointed several novel candidate nf-κb targets in these and other pathways that have not yet been described to be regulated by nf-κb. this is exemplified by identification of κb sites for diablo, which potentiates some forms of apoptosis, and for the traf family-associated nf-κb activator tank that inhibits nf-κb activation and enhances apoptosis by activating cell death signals and inhibiting survival signals. although it remains to be determined if these and other genes identified using ohmm are genuine nf-κb transcriptional targets, if confirmed, these could yield important new insights into the less well-characterized but nonetheless important pro-apoptotic activity of nf-κb that has been observed in response to certain stimuli and in certain cells. selected cellular pathways, biological functions and diseases in which our predicted nf-κb targets were over-represented are shown. the associated predicted nf-κb targets are represented by official human gene symbols. genes containing κb sites with predicted occupancy probability greater than were used in this analysis. please see additional file for the complete list. genes known in the literature to be regulated by nf-κb are denoted with *. similarly, ohmm pinpointed κb sites for ikbkb, which is the key nf-κb activating kinase in the canonical nf-κb signaling pathway, as well as for toll-like receptor that serves as a coreceptor for rna-associated autoantigens and participates in the innate immune response to microbial agents, potentially extending the list of nf-κb-regulated tlrs. if validated as genuine nf-κb targets, identification of ikkbeta and tlr by ohmm may have uncovered novel positive feedback loops for amplification of nf-κb signaling. our analysis suggests that nf-κb could play new roles in some of the pathways in which it is already known to participate, as seen in the protein ubiquitination pathway that can either regulate protein activation/function or target proteins for proteolytic degradation via the proteasome. for instance, while it is known that nf-κb activates the expression of deubiquitinating enzymes cyld and a that negatively regulate nf-κb signaling as well as expression of the proteasome subunit lmp and proteasome activators pa28-alpha and pa28-beta, ohmm identified the e ubiquitin-conjugating enzymes ube2h, ube2d and ube2m that promote protein ubiquitination, as putative nf-κb targets. incidentally, ube2m was previously observed in microarray studies to be induced by epstein-barr virus that leads to increased nf-κb activity, although it remains to be determined if ube2m is directly controlled by nf-κb. interestingly, our analysis also pinpointed κb sites in btrc. btrc is the substrate recognition component of a skp1-cul1-f- box protein e ubiquitin ligase complex that mediates ubiquitination and subsequent proteasomal degradation of target proteins including the nf-κb inhibitor proteins ikappab alpha, ikappab beta and ikappab epsilon, allowing free nf-κb dimers to translocate to the nucleus and activate transcription. btrc also mediates ubiquitination and subsequent proteasomal processing of the phosphorylated nfkb1/p and nfkb2/p precursor proteins that respectively participate in the canonical and non-canonical nf-κb signaling pathways. if btrc is a true target of nf-κb, this would suggest that upregulation of btrc by nf-κb could set up a positive feedback loop to amplify degradation of ikappab alpha, thereby providing a novel means for magnifying nf-κb signaling. overall, these results highlight a potential new role for nf-κb in promoting protein ubiquitination by modulating expression of e ubiquitin-conjugating enzymes and e ubiquitin ligases, thus expanding on its previously documented role in inducing expression of deubiquitinases, proteasome activators and proteasome subunits. ohmm also revealed new insights into the cross-talk between nf-κb and other signaling pathways. for example, we identified κb sites in several key mediators in the notch signaling pathway that is involved in cell-cell communications to regulate a broad spectrum of cell-fate determinations. these include the delta-like ligand for notch receptors dll the notch receptor, transcriptional regulator rbp that acts as a transcriptional repressor in absence of notch but is a transcriptional activator when associated with activated notch, and mastermind-like that serves as a transcriptional coactivator for notch. prior studies have suggested multiple cross-talk mechanisms between the notch and nf-κb signaling pathways, most of them suggesting activation of nf-κb downstream of notch, as well as protein-protein interactions between components of these two pathways. our finding that key mediators of notch signaling harbor κb sites supports the notion that nf-κb might in turn modulate notch signaling to influence cell fate determination during development, immunity and cancer. this would be consistent with recent work showing that nf-κb stimulates the expression of notch targets hes- and deltex- to stimulate marginal zone b-cell development and with prior observations from our group that nf-κb can trigger notch signaling by inducing expression of the notch ligand jagged although it remains to be determined if its effect on jagged expression is direct. other genes uncovered by ohmm could be intimately associated with nf-κb's role in transcriptional regulation, and in particular, with its recently uncovered role in gene-specific transcriptional repression. for example ohmm predicts κb sites for deacetylases hdac and sirt as well as for transcriptional corepressor sin3a. while nf-κb is known to engage in protein-protein interactions with histone deacetylases and corepressors to modulate gene transcription, ohmm raises the possibility of a new mode of action for nf-κb in this context. also highlighted were candidate targets in xenobiotic metabolism, including κb sites in sulfotransferase sult1c aldehyde dehydrogenase family gene aldh3b and transcription factor nuclear factor erythroid derived 2-like that regulates the oxidative stress response, in addition to the previously documented nf-κb target glutathione s-transferase. this points to a possible role for nf-κb in drug metabolism, multidrug resistance and detoxification of poisonous compounds, with possible impact for treatment of infectious diseases, anti-cancer therapy and/or environmental science. finally, ohmm pinpointed possible involvement of nf-κb in new pathways, such as the neurotrophin/trk signaling cascade, in which its role is emerging. in addition to the recently recognized nf-κb target phosphatidylinositol 3-kinase catalytic subunit pik3c, ohmm identified κb sites for its regulatory subunit pik3r as well as for the ras-related g-protein rras, ribosomal protein s kinase rps6ka and transcription factor cyclic amp-responsive element-binding protein. in this regard, creb was recently described to be controlled by nf-κb as is the critical neurotrophin nerve growth factor ngf that controls neuronal cell fate. these results may warrant further investigation in view of the emerging role of nf-κb in the regulation of neuronal survival and function in the nervous system, its effect on cognition and behavior and its suggested roles in epilepsy, stroke, alzheimer's, parkinson's diseases and other neurological disorders. verification of some of the predicted nf-κb sites in the human genome gel shift experiments were performed on several of the predicted nf-kb sites using two of the nf-κb family members: rela and c-rel. three low occupancy sequences were chosen as negative control and the known nfkb1a site was taken as the positive control. we chose the predicted target sites to be tested from some of the top categories that were extracted by pathway analysis. the sites to be tested were from the promoters of following genes: ikbkb and tlr involved in the nf-kb signaling pathway; diablo, tank, bcl2l and bnip involved in the death receptor signaling and apoptosis categories; hdac sirt and sin3a, involved in cell cycle g1/s and transcriptional regulation; and btrc, ube2d and ube2m, involved in ubiquitination pathway. all of these sites show stronger gel shift compared to the negative controls. to understand which whether predicted occupancy is a good indicator of binding compared to, say, evolutionary sequence conservation, we decided to classify these sites as the following as shown in figure the conservation score greater than 90% is considered high, between 90% and 80% is considered to be medium and less 80% is considered to be low. similarly, sites with occupancy greater than is considered to be high occupancy and the rest are considered to be low occupancy. the tested sites, organized according to this classification, are listed below. high conservation and high occupancy: ifbkb, ube2d and sin3a high conservation and low occupancy: tank- tank- bcl2l and ube2m medium conservation and high occupancy: sirt1 medium conservation and low occupancy: tlr hdac and btrc low conservation and high occupancy: bivip and diablo in general, the high occupancy sites are better binders, the exceptions being sin3a, hdac and tlr. for hdac the occupancy is only slightly lower than the cut-off threshold. in comparison, most highly conserved low occupancy sites have weak gel shifts, suggesting ohmm occupancy measure is a better predictor of binding affinity than the sequence conservation score. discussion while the use of an hmm is thus advantageous in identifying sites, we have discovered that an hmm offers tremendous benefits for the special case of identifying self-overlapping motifs. when putative sites overlap, the overall occupancy probability in a position in the dna sequence depends upon the strength of all sites containing that position. calculation of the overall occupancy probability over the entire sequence using the gamma variables in the hmm method enables us to compute the occupation probability in presence of overlapping sites. in this regard, our approach has a distinct advantage over the two hmm scoring methods in vogue at the present time, viz. the likelihood method and the viterbi method. the likelihood of a sequence scored over a long stretch of sequence is greater when more sites are present in that sequence. it is perfect for discovering clusters of sites, as evident from the work on fly regulatory modules. the likelihood score by itself does not indicate the location of sites in a promoter. in contrast, the viterbi method proclaims the presence of sites at positions where the state path with the highest probability contains the motif state. it, however, fails to consider overlapping sites because only the best site would be present in the optimal state path. training ohmm is a non-trivial two-step process. in the first step, we trained both the emission and transition probabilities using short sequences rich in known sites with the aim of accurately estimating the motif profile. in the second step, we kept the emission probabilities constant, and trained the transition probabilities on promoters containing known sites as well as random promoters to accurately estimate the transition probability to the motif z reflecting the site density in the promoters of all genes in the human genome. unlike most hmms in the literature, we successfully trained the hmm emission probabilities without requiring pre-alignment of training sites. one interesting feature is that trained z is inversely related to the strength of the motif profile. in other words, when the motif profile is kept constant and only z is trained, the weaker the motif profile used, the higher the trained z. this is probably the result of the compensating mechanism between the motif profile and z discussed above. interestingly, this competition between z and the motif profile also determines the occupancy probability as we saw before. the goal of the second training step was to estimate the transition probability to the motif z reflecting the site density in the promoters of all genes in the human genome. this z corresponds to the appropriate threshold when identifying sites in all human promoters. obviously, z trained in the first step was not appropriate due to the high site density in the site-rich training sequences. instead of training z on all human promoters, which is computationally expensive, we train it on a sufficient number of promoters to get a reasonable estimate of z. our procedure of beginning with the human promoters containing known sites as the training set and progressively adding promoters of randomly selected human genes to the training set leads to decreasing z, until the training set reached a few thousand promoters and z stabilized. we expect that, for a particular transcription factor, the speed of convergence of this procedure depends upon overall number of target genes of the factor: larger the number of targets/sites, quicker the convergence. one shortcoming of most site identification methods is that they assume equal probability of the site in a window of arbitrary length around the tss. the majority of known κb sites are located just upstream of the tss in gene promoters and the number of known κb sites decreases further upstream. specifically, of the known κb sites upstream of the tss, are located within bp and are located within bp of the tss. liu et al. have also made a similar observation in the promoters of nf-κb-regulated immune genes. to counter the claim that such an observation for binding sites may be due experimental bias, we cite tabach et al. who showed in a wide-scale bioinformatic study that functional binding sites are more likely to be present in the bp region upstream of the tss than any other upstream region for most human transcription factors and specifically for nf-κb. they defined functional binding sites as those over-represented in functionally related genes and conserved in related species. to bolster their conclusion, they showed location dependence of binding sites for the transcription factor myocardin in a controlled experiment. xie et al. also arrived at a similar conclusion based on binding site conservation. even though the exact reasons for the occurrence of such a phenomenon are not known at the present time, a better interaction of the transcription factor with the transcription machinery if it is bound close to the tss and the low density of nucleosomes near the tss heuristically explain why this phenomenon may occur. while we agree that a great number experiments need to be conducted to definitively prove the location dependence of functional κb sites, we feel that a site identification method must be able to take it into account. noticing that site density decreases sharply with the upstream disease from the tss and that transition probability to the motif z is proportional to site density, we modelled z using an exponential functional form such that the probability of identifying a site decreases further upstream but never reaches zero. even though this approach may fail to identify sites in distal promoters and enhancers, we believe that it allows site search in large upstream regions without identifying too many false positives. a κb site exerts influence on the occupancy probabilities at the positions surrounding the site in either direction. because overlapping binding sites are usually present when the scoring window is moved by one position, occupancy probability stays high at all positions in that window. when the window is moved by more positions, the occupancy probability at the new positions dips slightly below the average due to the high motif probability at the site. the occupancy probability returns to the average background value when the scoring window moves by ten positions. the relationship between the occupancy probabilities of sites located close to each other is quite instructive. when two sites are in tandem without any space between them, occupancy probabilities of both of them are lower due to the small motif-to-motif transition probability. occupancy probabilities are not very high even when the sites are one position apart because a window shift by one position from a κb site usually contains an overlapping κb site. the sites need to be at least two positions apart so that they do not exert significant influence on each other's occupancy probabilities. please note that in any case, the overall occupancy in a region containing two nearby sites is quite high. a strong correlation exists between the occupancy probabilities predicted by ohmm and the in vitro binding affinity of the nf-κb family members c-rel and rela for oligonucleotides in a gel shift experiment. this validates our physical binding model because the hmm-predicted occupancy probabilities appear to correspond to observed binding affinities. in the process of the study, we found how the padding sequence affects the results. the current practice for performing gel shift experiments in the nf-κb community consists of using particular padding sequences around the 10-mer corresponding to a potential κb site. the commonly used padding sequences, however, can form spurious κb sites due to the self-overlapping nature of κb motif, and hence the experiment may not represent binding of the 10-mer in the native promoter. ohmm offers guidance on the selection of correct padding sequences when designing experiments. a gel shift oligonucleotide should have an a 5' to the κb site and a t 3' to the κb site for minimum interference, and that a 3' c should be avoided at all cost. ideally, oligonucleotides containing a few bases corresponding to those surrounding the κb site in the promoter of the gene should be used to pad the site. this will capture the effects of all the overlapping putative binding sites in the native promoter and will avoid creation of artefacts based on the nucleotides present in the padding sequence. the nf-κb family of transcription factors consists of some of the most important proteins for an organism's development and survival. in resting cells, inactive nf-κb dimers reside in the cytoplasm, bound to inhibitors in the iκb family. many different stimuli can activate latent nf-κb/iκb complexes, including pro-inflammatory cytokines that trigger rapid activation of the classical nf-κb signaling cascade via the ikk kinase complex. consequent phosphorylation and degradation of iκb via the proteasome culminates in nuclear translocation of active nf-κb dimers and their binding to κb sites. this commonly results in transcriptional activation, or occasionally in repression, of genes important for the immune and inflammatory response, cell proliferation, adhesion, angiogenesis and inhibition of apoptosis. the non-canonical nf-κb signaling cascade is activated via the lymphotoxin and baff receptors, leading to activation of ikkα, nuclear translocation of p52/relb complexes and plays a role in lymphoid organogenesis and b-cell function. hundreds of nf-κb-regulated genes have been identified. among them, binding of activated nf-κb to the iκbα promoter can trigger its resynthesis, giving rise to an autoregulatory loop that terminates the activation process. consequently, under normal conditions, rel/nf-κb activation is tightly regulated and transient. however sustained activation of rel/nf-κb can promote cell survival, proliferation and/or cell death depending upon the nature and the extent of the stimulus, and is implicated in many pathological conditions including immune system disorders, chronic inflammation and cancer. although several nf-κb target genes have been identified, whether other nf-κb-regulated genes are its direct or indirect targets is not known. in addition to identifying many of these known nf-κb target genes, ohmm pinpointed several novel candidate nf-κb targets that are yet to be shown as controlled by nf-κb. while further studies will be needed to determine which of the candidate genes identified by ohmm are genuine transcriptional targets of nf-κb, these are likely to yield important new insights into the less well-characterized roles of nf-κb in the induction of apoptotic cell death, its recently uncovered role in transcriptional repression, the regulation of protein ubiquitination and its role in other signaling pathways, to name a few. in addition to identifying multiple genes with key roles in haematological and immune cell development and function, ohmm pinpointed κb sites in genes with roles in drug metabolism, skeletal and muscle system development and function, as well as in disease conditions ranging from immune and inflammatory disorders to infectious diseases, cancer, skeletal muscular disorders and neurological diseases. these findings are consistent with experimental evidence suggesting a role for nf-κb in immunological and inflammatory diseases, cancer and therapy-resistance, in skeletal myogenesis and cachexia, as well as in cognition, behavior and neurological disorders like epilepsy, stroke, alzheimer's, parkinson's diseases, and amyotrophic lateral sclerosis. ohmm is thus a unique tool to shed new light on the multifaceted biological functions of nf-κb and the mechanisms involved, and to identify critical mediators of dysregulated nf-κb activity in human disorders that may ultimately help to develop new approaches for treatment. let us finish the discussion by mentioning possible shortcomings of the method. ohmm requires a complicated two-step training procedure, which might only work practically for factors with many targets. we have not considered insertions or deletions inside the motif, making ohmm well suited for the simultaneous training of emission and transition probabilities. in addition, we have considered the transition probability to the motif z in the upstream region to decrease exponentially during the scoring procedure. this density function may not be accurate, especially for enhancers. we will need many more sites to establish a more accurate density function.  our method is the first attempt to deal specifically with self-overlapping binding motifs, providing guidance in the selection of padding sequences in gel shift experiments. when the hmm is considered as a biophysical model, it interprets the weight matrix as binding energy, the transition probability to the motif as transcription factor concentration and the gamma variable as occupancy probability. because it uses occupancy probability as the discriminant function, it learns the associated natural threshold in a principled manner during the training procedure. another unique feature of ohmm is that transition probabilities change with sequence location to reflect site density. transcription factor's occupancy probability on a site can be accurately estimated using our method, as seen by the high correlation with experimental binding affinities. high evolutionary conservation scores and enrichment in experimentally regulated genes suggest that κb sites predicted by our method might be functional. our results may provide important new insights into the function and regulation of nf-κb and uncover possible new biological roles for this important transcription factor family. 
while identifying sites, it adjusts parameters to model site density changing with the distance from the transcription start site. ohmm could be used to predict physical occupancy, and provides guidance for proper design of gel-shift experiments. ignoring these features can lead to inaccurate predictions as well as incorrect design and interpretation of experimental results. moreover, considering ohmm as a biophysical model allows us to learn the classification threshold in a principled manner. high evolutionary conservation scores and enrichment in experimentally verified regulated genes suggest that nf-κb binding sites predicted by our method are likely to be functional. based upon our predictions, new insights into nf-κb function and regulation and possible new biological roles of nf-κb were uncovered. interpreting hmm parameters as biophysical quantities, this method uses the occupancy probability of a transcription factor on a dna sequence as the discriminant function, earning the algorithm the name ohmm: occupancy via hidden markov model. in the context of binding sites to transcription factor nf-κb, we find that the occupancy probability predicted by ohmm correlates well with the binding affinity in gel shift experiments. dna sequence binding motifs for several important transcription factors happen to be self-overlapping. we have developed a tool based on a hidden markov model that identifies binding location of transcription factors with preference for self-overlapping dna motifs by combining the effects of their alternative binding modes. our method deals specifically with identifying locations with multiple overlapping binding sites by computing the local occupancy of the transcription factor. ohmm learns the classification threshold by training emission probabilities using unaligned sequences containing known sites and estimating transition probabilities to reflect site density in all promoters in a genome. another feature of ohmm is that we allow transition probabilities to change with location relative to the tss. moreover, most methods use arbitrary thresholds and fail to provide a biophysical interpretation of statistical quantities.
 dna sequence binding motifs for several important transcription factors happen to be self-overlapping. many of the current regulatory site identification methods do not explicitly take into account the overlapping sites. moreover, most methods use arbitrary thresholds and fail to provide a biophysical interpretation of statistical quantities. in addition, commonly used approaches do not include the location of a site with respect to the transcription start site in an integrated probabilistic framework while identifying sites. ignoring these features can lead to inaccurate predictions as well as incorrect design and interpretation of experimental results.  we have developed a tool based on a hidden markov model that identifies binding location of transcription factors with preference for self-overlapping dna motifs by combining the effects of their alternative binding modes. interpreting hmm parameters as biophysical quantities, this method uses the occupancy probability of a transcription factor on a dna sequence as the discriminant function, earning the algorithm the name ohmm: occupancy via hidden markov model. ohmm learns the classification threshold by training emission probabilities using unaligned sequences containing known sites and estimating transition probabilities to reflect site density in all promoters in a genome. while identifying sites, it adjusts parameters to model site density changing with the distance from the transcription start site. moreover, it provides guidance for designing padding sequences in gel shift experiments. in the context of binding sites to transcription factor nf-κb, we find that the occupancy probability predicted by ohmm correlates well with the binding affinity in gel shift experiments. high evolutionary conservation scores and enrichment in experimentally verified regulated genes suggest that nf-κb binding sites predicted by our method are likely to be functional.  our method deals specifically with identifying locations with multiple overlapping binding sites by computing the local occupancy of the transcription factor. moreover, considering ohmm as a biophysical model allows us to learn the classification threshold in a principled manner. another feature of ohmm is that we allow transition probabilities to change with location relative to the tss. ohmm could be used to predict physical occupancy, and provides guidance for proper design of gel-shift experiments. based upon our predictions, new insights into nf-κb function and regulation and possible new biological roles of nf-κb were uncovered. 