BACKGROUND
identifying mechanisms that regulate gene expression in bacteria is essential for understanding and eventually controlling their pathogenicity.
all known bacteria share a well conserved transcriptional holoenzyme, rna polymerase .
the rnap is comprised of a 3-subunit catalytic core plus a variable σ-factor subunit that provides dna binding specificity.
one of these σ-factors, σ <dig> in escherichia coli, participates in the transcription of a majority of genes including those with housekeeping functions.
e. coli is the best studied bacterial model with regard to promoter identification and prediction.
as such, most promoter predictions are based upon the analysis of e. coli σ <dig> promoter data.
the earliest collections of e. coli σ <dig> promoters revealed the - <dig> and - <dig> hexamer consensus motifs, ttgaca and tataat, that serve as recognition sites for the  <dig>  and  <dig>  domains of σ <dig>  <cit> .
position weight matrices  were the first models to quantify the hexamer motifs  <cit> .
pwm models were expanded to quantify the variable-length spacer region between hexamers  <cit> , which is important in orienting the hexameric motifs for interaction with the sigma binding factors  <cit> .
challenges encountered by pwm models include defining thresholds that are sensitive enough to include known promoters without predicting numerous false positives.
most of the quantitative modeling efforts that ensued require training sets comprised of both positive and negative sequences.
artificial neural networks   <cit>  have been trained on sequences of identified e. coli promoters and non-promoters.
a hidden layer in the ann architecture quantifies interactions among pairs and triplets of nucleotides.
the resulting ann scans and scores overlapping sequences, and reports a score in the range  that indicates the likelihood of the sequence being a promoter.
a time-delay neural network  can combine two simple anns  with a variable-length spacer region  <cit> .
burden et al   <cit>  measured the distance from the transcription start site  to the translation start site  of  <dig> e. coli promoters.
they showed that the distribution peaks sharply around  <dig> nt, and that combining the tss-tls distribution with the nnpp <dig>  tdnn  <cit>  significantly enhances the specificity of the prediction.
in another machine learning approach that has been applied to model promoters, support vector machines  were trained on e. coli promoter sequences of length  <dig>  <cit> .
although the svm approach has the advantage of comprehensively quantifying the primary structure of the upstream region, it does not examine structures of higher order that motivate our approach.
a natural extension to pwms that explicitly models an empirical spacing distribution between motifs is given by duration hidden markov models .
here "duration" refers to this explicit representation of a spacer length distribution, as opposed to the geometrically distributed lengths that are expected from components of profile hidden markov models  <cit> .
although the variable-length spacer region between hexamers has been incorporated into promoter modeling and predictors before  <cit> , none of these earlier efforts have integrated an explicit probabilistic representation of the spacer distribution into a reusable predictor as a duration hmm, which is arguably its most natural representation.
on the other hand, while duration hmms have been introduced into genome analysis , they have not to our knowledge been applied to modeling transcriptional or translational signals before.
bacteria of the genus chlamydia are obligate intracellular parasites that were genetically isolated from other bacteria nearly a billion years ago when they moved into their intracellular environment  <cit> .
in humans, chlamydia infections are responsible for infertility, blindness, arthritis and cardiovascular disease  <cit> .
because chlamydiae have an intracellular life-cycle, standard genetic techniques are often insufficient to study gene regulation  <cit> .
hence, only  <dig> to  <dig> promoters have been experimentally verified  <cit> .
however, with a small genome of only about  <dig> mbp and  <dig> genes, chlamydia trachomatis  makes a good candidate for in silico analysis.
surveys of known bacterial promoters suggest that their structures are relatively diverse  <cit> .
additionally, established ct promoters display obvious differences from the established consensus hexamers of e. coli  <cit> .
although σ <dig>  the ct analog of e. coli σ <dig>  has dna binding domains homologous to domains  <dig>  and  <dig>  in σ <dig>  sequence based phylogenetic analysis of bacterial rnap subunits has shown discernable evolutionary distance between the ct and e. coli rnap subunits  <cit> .
therefore, it is plausible that an organism-specific model is appropriate for ct.
phylogenetic footprinting takes advantage of relative conservation of motifs among related species.
grech et al   <cit>  developed an algorithm that combines e. coli trained pwms and chlamydial phylogenetic footprinting.
ct upstream regions are screened with the pwms and the potential promoter hexamers are filtered with an algorithm that accepts only conserved sequences in a consensus of c. trachomatis, c. pneumoniae and c. caviae.
although this is a promising approach, because they used an e. coli trained pwm, their results may be strongly influenced by prior expectations that all bacterial promoters are structured as in e. coli.
we believe that more development is needed in ab initio approaches for predicting promoters using sequence information directly from the organism under study  in combination with biophysical metrics that derive from known models about the biology of transcription in general.
this study aims to develop ct promoter models using only known ct promoters.
to do so, it considers dna stability and topological features of the upstream region as well as rnap σ-factor/dna binding.
as hertz and stormo   <cit>  aptly wrote "... the polymerase needs to bind the dna, open the dna, initiate transcription, and release the promoter for elongation."
the tdnns and svms that consider extended promoter sequences are addressing this issue from a sequence perspective.
this study utilizes measures that have been developed to quantify stability and other aspects of dna structure.
evidence from the profiling of dna curvature, bendability, twist, stability and propensity for stress-induced destabilization in e. coli, b. subtilis, c. trachomatis, plants and vertebrates  <cit>  suggests that there are peaks for these measures near the tss.
here we use a stochastic model building procedure that allows for the combination of relevant predictive measures selected from rnap σ-factor/dna binding propensity, as quantified by duration hmms, and structural features of the upstream region, as quantified by biophysical metrics.
methods
stochastic model building
stepwise binary logistic regression   <cit> , as implemented in spss version  <dig>  statistical software , selects an optimal set of independent variables  to classify observations into two populations.
logistic regression does not assume a linear relationship between the dependent and independent variables, normal distributions, or homoscedasticity .
it does, however, assume independence of observations.
we address this requirement in a separate section describing the selection of non-redundant observations.
the mathematical model  fitted by sblr has the form
where i is the number of steps, v <dig> through vi are the predictor variables selected, and b <dig> through bi are coefficients determined by the analysis.
u is the logit for the dependent variable, which means that
here, the event is class membership.
when p denotes the prob, the equation can be rewritten as
selecting a cutoff for p, most commonly  <dig> , converts p into a classifier.
when  <dig>  is the probability threshold, eu =  <dig> and the classification threshold for u is  <dig>  the effectiveness of a model can be evaluated by its ability to correctly classify the training data.
the spss sblr analysis procedure provides many user-defined options.
we selected the forward conditional stepwise procedure for all analyses.
at each step, a score statistic is calculated for each variable excluded from the model.
the score statistic is based on maximum-likelihood estimation criteria and is asymptotically distributed as a χ <dig> variable  <cit> .
the variable with the highest significant χ <dig> value is entered into the model.
if no significant variables remain, then the procedure stops with the current model.
similarly there is a mechanism for stepwise removal.
after a new model has been generated, score statistics are calculated for all variables in the model.
if the p-value for any variable in the model is greater than the probability for stepwise removal, then the variable is removed from the model.
we retained the default probabilities for stepwise entry  and removal , thus ensuring that the significance of all model variables is less than  <dig> .
potential observations and dependent variable
a significant challenge for bioinformaticians is to model data that has been collected by multiple laboratories using different assays, protocols and equipment.
this phenomenon is compounded in the study of ct where the organism is metabolically active only inside an infected host-cell.
one way to minimize the use of conflicting and/or controversial data is to rely upon reviews written by informed biologists.
for this reason, we consulted the reviews of mathews & timms   <cit>  and tan   <cit>  to compile a list of  <dig> experimentally verified σ <dig> promoters.
to this list we added  <dig> promoters that were experimentally verified by grech et al   <cit>  and hefty et al   <cit>  after the previously cited reviews were written.
for the purposes of this study, we consider these  <dig> sequences to be the known ct σ <dig> promoters.
ant distance from the lead nt of the - <dig> hexamer to the tls.
breferences: m: mathews & timms  <cit> ; g: grech et al  <cit> ; h: hefty et al  <cit> ; t: tan  <cit>
chour post infection of transcriptional activation  <cit>
<dig>  files containing the ct genome  and genome table  were retrieved from the ncbi website, , on july  <dig>   <dig> and last modified by ncbi on january  <dig>   <dig>  an r script was written to extract the  <dig> nt upstream regions of all  <dig> protein-coding genes as annotated in the genome table.
the regions were verified at several stages throughout the research with other sources, including a ct genome database previously available from los alamos national laboratories and comparable prediction algorithms.
table  <dig> displays the distance from promoter to tls for all training set genes.
since the maximum distance is  <dig> nt,  <dig> nt was set as the upper limit for data analysis.
then, the upstream region was defined as  <dig> nt to allow for biophysical structures  <dig> nt upstream from a predicted promoter.
<dig>  for each of the  <dig> genes listed in table  <dig>  the  <dig> nt upstream region was parsed into overlapping sliding window sequences of length  <dig>  and step-size  <dig>  each subsequence  was labeled according to its parent gene and position occupied in the upstream region: e.g.
the first seq <dig> was labeled ct046_ <dig> because the initial nt is found  <dig> nts upstream from the ct <dig> tls.
<dig>  the dependent variable, promoter, was assigned a  <dig> if a promoter sequence listed in table  <dig> was totally contained in seq <dig>  and  <dig> otherwise.
thus, 1's identify potential promoter observations and 0's identify potential non-promoter observations.
<dig>  cases with upstream positions ≥  <dig> and ≤  <dig> were selected as potential observations to restrict the analysis to the range of the training set data.
the upper bound is  <dig> nt upstream of the furthest upstream training set promoter and the lower bound is equal to the furthest downstream training set promoter.
independent variables
the primary variable for promoter prediction is the pattern that characterizes the binding between the rnap σ-factor and dna.
here we use duration hmms to describe and quantify rnap-σ66/dna binding.
after a set of known promoters is used to train a duration hmm, the duration hmm scans a new sequence to identify the hexamer-spacer-hexamer subsequence that scores the highest with regard to potential binding.
the variable hmm_score is assigned the score associated with the highest scoring subsequence, while the variable start denotes the position of the lead nucleotide in the - <dig> hexamer and end denotes the position of the last nucleotide of the - <dig> hexamer.
specifically, a training set of promoter sequences was placed in the file ts.txt.
the initial ts.txt contained the contents of table  <dig>  columns "- <dig> hex", spacer and "- <dig> hex".
this file was supplied as input to durahmmer  which was used to create a duration hmm with the command: durahmmer - <dig>  <dig> - <dig>  <dig> -s  <dig> -s  <dig> -p  <dig> -u  <dig> : <dig> : <dig> : <dig>  -c ts.txt > ts.hmm.
the options to the command specify the following model parameters:  <dig> matched states  at the 5' and 3' sequence ends; minimum and maximum spacer lengths of  <dig> and  <dig> respectively; a background compositional model of  <dig> % a,  <dig> % c,  <dig> % g, and  <dig> % t; and spacers should be modeled to have their empirical composition in the training set .
the program durahmmer produces a valid hmmer  <dig> .
<dig>  <cit>  model file representing a duration hmm.
for the final model of this study, the model file and the input data file are provided as additional files  <dig> and  <dig>  all  <dig>  seq <dig> observations from the  <dig> genes were placed in the file all.txt so that optimal promoters and hmm scores could be calculated by hmmsearch  <cit>  with the command: hmmsearch -e  <dig> ts.hmm all.txt.
we ran hmmsearch with a high e-value because we were interested in combining the score of the maximum scoring hit with other metrics in a composite procedure regardless of its magnitude.
in combination with the duration hmm model score described above, we also used the following biophysical metrics of promoter position and structure as possible independent variables for the sblr model:
<dig>  position, which indicates the location of seq <dig> in the upstream region relative to the tls.
for ct046_ <dig>  position =  <dig>
<dig>  measures of curvature   <cit>  and %gc content  for each  <dig> nt upstream region, which were determined by the online bend.it server  with a window-size of  <dig>
<dig>  free energy change  of dna melting , bendability  and twist angle , which were determined for each  <dig> nt upstream region by the online plot.it server .
all measurements were then averaged over each seq <dig>  Δg always has a negative sign and is interpreted as greater values having lower stability.
for statistical analysis this variable was transformed by stable = -Δg so that the sign is always positive and the interpretation is that larger values have greater stability.
stability is also of interest in the immediate downstream region, so positions 27- <dig>  and 1- <dig>  were quantified.
since the bendability measure increases with rigidity, it was renamed rigid.
the twist angle measurement, twist, was not transformed.
<dig>  possible times of expression onset include  <dig>   <dig>   <dig>   <dig> and  <dig> hours post infection .
mutually exclusive binary variables h <dig>  h <dig>  h <dig>  h <dig>  h <dig> and h <dig> were created to mark time of expression onset.
<dig>  stress-induced dna duplex destabilization  quantification utilizes structural and energetic properties of dna to measure the propensity for strand separation under negative superhelical stress  <cit> .
a low sidd score indicates a high propensity for strand separation.
sidd measurements were determined by the websidd server  <cit>  with default parameters except for open region size =  <dig>  because niehaus et al  <cit>  have shown a time dependent response to chlamydial dna supercoiling, interactions between the time of expression onset and sidd were included  <cit> .
the sidd/hour of onset interaction is quantified by sidd_h# = sidd*h#.
<dig>  for variables based on the entire seq <dig>  lagged variables were created for the four non-overlapping upstream subsequences of length 32: e.g.
for ct046_ <dig>  curve_l <dig> was set equal to the curve value of ct046_ <dig>  curve_l <dig> was set equal to the curve value of ct046_164; curve_l <dig> was set equal to the curve value of ct046_196; and curve_l <dig> was set equal to the curve value of ct046_ <dig>
selection of non-redundant observations from potential observations
as mentioned earlier, sblr assumes independent observations.
to address this requirement, we select for analysis a subset of the overlapping potential observations that are non-redundant with respect to the pair of hexamers that are most likely to bind the rnap σ-factor.
if we select only those cases where end =  <dig>  we eliminate all of the redundant optimal hmm hexamer pairs while retaining most optimal hmm instances .
table  <dig> demonstrates how this selection ensures that neighboring optimal hmm instances that match are included only once.
six potential observations, ct046_ <dig> through ct046_ <dig>  all contain the verified promoter with hexamer pair tggtta and taaaaa.
consequently, they all have promoter =  <dig> and hmm_score = - <dig> .
but only ct046_ <dig> has end =  <dig> and is selected to represent the verified ct <dig> promoter.
similarly, only ct046_ <dig> represents the maximal non-promoter hexamer pair ttgtgt and aaaagt with score = - <dig> .
this process incidentally aligns each selected seq <dig> such that the optimal downstream hexamer is at the far right end.
this selection process does not eliminate overlapping sequences, but it does eliminate overlapping likely binding sites.
ct046_ <dig> and ct046_ <dig> overlap a great deal.
however, the last hexamer of ct_046_ <dig>  is not present in ct046_ <dig> and the first hexamer of ct_046_ <dig>  does not appear in ct046_ <dig>
it should be noted that although each training set gene begins with the same number of potential observations, this selection process causes the number of selected non-redundant observations to differ among genes.
each gene starts with around  <dig> potential observations with promoter =  <dig> for each verified promoter, and around 325-40- <dig> =  <dig> potential observations with promoter =  <dig>  however, selection for non-redundant observations always results in the number of designated non-promoters being reduced to approximately  <dig>
while selecting sequences with non-redundant hmm_scores does mitigate the problem of dependent observations, it may not entirely eliminate it.
while there are numerous studies that affirm the robustness of baysian discriminant analysis with regard to violating the assumptions of a linear relationship between the dependent and independent variables, normal distributions, and homoscedasticity  <cit> , we could not find similar studies regarding the robustness of logistic regression.
an alternative to the current analysis would be to use stepwise discriminant analysis, knowing that we are violating some assumptions.
there are versions of logistic regression, including generalized estimating equations   <cit> , that are specifically designed for correlated data such as longitudinal studies.
in these procedures there are subject variables and within-subject variables.
it might be possible to force this study data into such a format, but as yet there are no readily available stepwise procedures to scan multiple possible predictors.
a final alternative would be to select non-overlapping sequences with the penalty of losing information and perhaps introducing a selection bias.
sblr is a procedure for model identification.
it is only after a model has been identified that it can be evaluated for independence.
given that, we elected to analyze the non-redundant observations with sblr and then examine the error terms for independence.
in time series analysis , this is done by checking that the error term is normally distributed with zero mean, and that autocorrelations and partial autocorrelations of the error term are not significant  <cit> .
iterative modeling strategy
sources of error that could lead to misclassification include  imprecise laboratory procedures in defining and identifying promoters ,  presence of more than one promoter population,  failure to include relevant predictor variables, and  random variation.
to minimize the first two error sources, an iterative strategy was developed.
duration hmm iteration  addresses error source , while sblr iteration  addresses source .
duration hmm iteration
minor modifications in the configuration of the training set promoters can improve classification accuracy.
to accomplish this, we allowed each promoter to vary within a neighborhood that extends the sequence by  <dig> nts on each side.
a limit of  <dig> nts ensures that a modified hexamer will not locate completely outside of the original promoter sequence.
for example, when the promoter ct <dig> is extended, it becomes ttgtttgcagagtttttattttaaatatgttataatctgtc, with the bolded nts marking the extensions.
initially, a duration hmm is determined by the original, non-extended, promoter training set.
then the set of extended promoters is searched for the highest scoring instance of the duration hmm in each extended sequence.
if a high-scorer is not the same as the original promoter, it replaces the original in the training set.
the iteration continues until stabilization.
for the final model, ct <dig> was modified to ttgcagagtttttattttaaatatgttataat.
sblr iteration
deletion and subsequent replacement of members of the training set can eliminate promoters that are likely to be members of a different promoter population.
this is accomplished via the iterative scheme diagrammed in figure  <dig>  initially, the complete set of  <dig> verified promoters determines the duration hmm and the independent observations selected for sblr analysis.
sblr delivers a mathematical model that produces a predicted probability of class membership  for each observation.
a threshold on p of .
<dig> is used to classify each observation as a predicted promoter or non-promoter.
for those  <dig> cases where promoter =  <dig>  we also use the value of p to determine when a promoter appears to be an outlier and should be eliminated from the training set.
after observing the  <dig> probabilities, a retention threshold on p between  <dig> and .
<dig> is established.
if a training gene has only one identified promoter and that promoter has a p less than the retention threshold, then all observations for that gene are deleted from the analysis.
similarly, if a training set gene has two identified promoters and they are both selected for deletion, all observations for that gene are deleted.
however, if a training set gene has two identified promoters and only one is selected for deletion, all upstream observations for that gene remain in the analysis dataset and only observations within the remaining promoter are assigned promoter =  <dig>
modifying the training set in any way necessitates the determination of a new duration hmm, which in turn determines which observations will be aligned such that end =  <dig> and subsequently included in the next sblr analysis.
the iteration process continues until the training set stabilizes.
stratified k-fold cross-validation
once the final training set and model are selected, it is necessary to validate the model to ensure against over-fitting and to allow for comparisons with algorithms trained on other datasets.
in the case of dichotomous classification, stratified k-fold cross-validation  <cit>  partitions the training set into k subsamples such that each subsample has approximately the same proportions of class membership.
here we designate each training gene as a subsample; hence k equals the number of genes in the training set.
then, one gene  is retained as a validation set while the remaining genes are used as training data.
evaluation measures are calculated by aggregating the results of each validation set.
comparable algorithms
the following three algorithms were used to compare performance and to identify co-predictions with the model developed in this study: nnpp <dig> , tss-predict, and footy.
nnpp <dig>   <cit>  is an online time-delay neural network that is accessible for promoter predictions at .
we used the following options: organism = prokaryote and minimum promoter score =  <dig>  to define promoters in the  <dig> nt upstream region of all ct genes.
for the support vector machine algorithm tss-predict  <cit> , the top two ranking predictions for each ct gene are posted as supplementary material at doi: <dig> /j.combiolchem.
<dig> .
<dig>  the  <dig> ct promoters predicted by footy  <cit> , an algorithm that utilizes phylogenetic footprinting, are reported directly in the publication that describes the algorithm.
r scripts scanned the promoters predicted by nnpp <dig>  and tss-predict for matches with the promoters predicted by the study model.
an nnpp <dig>  match was declared when the study prediction was contained within the  <dig> nt nnpp <dig>  prediction.
a tss_predict match was declared when the tss_predict predicted hexamer pair was contained within the study prediction.
RESULTS
finding the best model
the initial model, m <dig>  utilizes the initial training set of  <dig> promoters with observations from their  <dig> parent genes.
the duration hmm model converged after one iteration, modifying the alignment of  <dig> promoters.
for all models, table  <dig> reports the variables that were selected for the model and evaluation measures.
if tp = true positive, fp = false positive, tn = true negative and fn = false negative, then sensitivity or recall = tp/, specificity = tn/, positive predictive value  or precision = tp/, negative predictive value  = tn/, and accuracy = /.
the total number of observations for each model differs according to the promoter training set being used.
training set
ct681a
ct681b
ct681a
ct681a
+stable1_37
-position
+curve_l32
-gc_l128
+rigid_l96
+stable1_37
-gc_l32
-position
+curve_l32
-curve_l64
-gc_l128
+stable1_37
-position
+curve_l32
-stable_l32
+sidd_h24
-curve_l128
-sidd_l128
+stable1_37
-position
+curve_l32
-stable_l32
-stable27_37
ppv or
athe variables are listed in order of entrance into the model and the sign indicates the sign of the coefficient.
broc analysis area under the curve
for model m <dig>   <dig> of the  <dig> promoters were classified correctly, with  <dig> false positives.
there is always the possibility that these are yet to be recognized promoters, but at this point they are counted as misclassifications.
for the  <dig> verified promoters that were missed, the predicted probabilities ranged from  <dig>  to  <dig> .
since a natural separation appeared to between  <dig>  and  <dig> , p =  <dig>  was selected as the retention threshold and promoters ct <dig>  ct681a, ct681b and ct <dig>  were deleted from the training set for the next model, m <dig>
the duration hmm model for m <dig> converged after two iterations, modifying the alignment of  <dig> promoters.
table  <dig> shows that m <dig> classified the modified training set perfectly, indicating that perhaps too many promoters had been deleted from the original training set.
the retention threshold was reset to  <dig>  and ct <dig> was reinstated for model m <dig>
the duration hmm model for m <dig> converged after one iteration.
table  <dig> displays the alignments of the  <dig> promoters that were modified.
m <dig> also classified the modified training set perfectly.
again the results indicated that the next model, m <dig>  should reset the retention threshold to  <dig>  and reinstate ct <dig>  however, table  <dig> reports that m <dig> is not as good as models m <dig> and m <dig> because of classification errors.
given two models, one training set a subset of the other, that both classify their respective training sets with 100% accuracy, we reasoned that the model trained on the largest set would provide the most sensitive genome-wide prediction.
thus, m <dig> was selected as the best and final model because of the perfect classification with the largest training set.
the complete data file used to build m <dig>  additional file  <dig>  is supplied so that others may replicate or modify the model.
finally, the error terms of m <dig> were checked for independence.
residuals, promoter - p, were calculated for all selected observations and shown to be normally distributed with zero mean.
additionally, the autocorrelations and partial autocorrelations of the residuals were not significant.
thus, the independence assumption of sblr was not violated by this model.
aggregated results of the stratified k-fold  m <dig> cross-validation are reported in the last column of table  <dig>  for the  <dig> genes in the m <dig> training set,  <dig> promoters  were not identified  and there were  <dig> false-positive predictions .
the incorrect classifications are most likely due to incomplete representation of the sample space, but may indicate additional populations or absent predictors.
model interpretation
the m <dig> duration hmm describes and quantifies the rnap-σ66/dna binding observed in the training set.
a visualization of the m <dig> parameters is shown in figure  <dig>  the - <dig> hexamer is dominated by the initial ttg motif, while the initial t with frequent as and ts describe the - <dig> hexamer.
the c and g compositions  of the spacer region are much smaller than those of the genome .
spacer lengths of  <dig> predominate, while spacers of length  <dig> are absent.
the input data file for durahmmer  and the resulting output data file  are provided as additional files  <dig> and  <dig>  the output data file is an hmmer  <dig> .
<dig> model file which supplies the parameters of the m <dig> duration hmm to hmmsearch.
complete documentation for the contents of the file can be found in the hmmer user's guide at .
briefly, the first  <dig> lines are header information with the main model section following.
there are  <dig> lines for each of the  <dig> possible nodes.
the first and last  <dig> nodes refer to the - <dig> and - <dig> hexamers, while nodes  <dig> through  <dig> refer to possible spacer positions.
the first line for each node displays the contribution to the final score  for the corresponding nucleotide matching a, c, g or t. the third line is particularly relevant to nodes  <dig> through  <dig>  which correspond to spacer nucleotides  <dig> through  <dig>  as nucleotides in these positions may or may not be present in the sequence being scored due to variable spacer length, the third line provides the odds of transitioning to another spacer nucleotide or to the - <dig> hexamer.
the m <dig> prediction equation generated by sblr is:
being the strongest predictor, hmm_score is selected in the first step of the sblr procedure.
the prediction equation for step one is
using a classification cutoff of p =  <dig>  and setting u =  <dig> yields hmm_score =  <dig>  as the threshold for step  <dig> classification.
at step  <dig>  14/ <dig> promoters and 2220/ <dig> non-promoters were classified correctly.
thus, the remaining eight model variables moved  <dig> promoters with hmm_score < .
<dig> to promoter classification and  <dig> non-promoters with hmm_score ≥ .
<dig> to non-promoter classification .
the predictor variables and their coefficients describe the verified promoters and their upstream regions.
promoters have high hmm_score and low position.
the near upstream region is curved and unstable, whereas the further upstream region is uncurved and unstable under superhelical stress.
for late-cycle genes where expression onset occurs at  <dig> h pi, the effect of superhelical stress is less than at other times .
the upstream characteristics may reflect transcription factor binding and/or additional interaction with the rnap holoenzyme.
the interpretation of the positive coefficient for stable1_ <dig> is more subtle.
in the second step of the sblr, four observations change from fp to tn and  <dig> observations change from fn to tp.
the means of stable, stable1_ <dig> and stable33_ <dig> are all larger in the second group than in the first.
although stable33_ <dig> shows the greatest mean difference, the most statistically significant is stable1_ <dig>
model exercise: predicting promoters for the ct genome
finally, the m <dig> model was used to predict promoters for the entire ct genome.
additional file  <dig> reports  <dig> predicted promoters in  <dig> unique genes, along with their hmm scores and genome locations.
thus, for  <dig> of the total  <dig> ct genes, this model does not find any 32-mers with a probability >  <dig> .
this suggests a conservative prediction that emphasizes specificity over sensitivity.
other explanatory factors may include alternate binding patterns for σ <dig>  alternative σ-factors, and operon configurations.
there was a substantial overlap among predictions by different methods.
additional file  <dig> lists the  <dig> promoters  co-predicted by m <dig> and nnpp <dig> , while additional file  <dig> lists the  <dig> promoters  co-predicted by m <dig> and tss-predict.
additional file  <dig> reports the  <dig> promoters  co-predicted by m <dig>  nnpp <dig>  and tss-predict.
all predictions are for  <dig> = position =  <dig>  consistent with the range of the modeling procedure.
of the  <dig> promoters predicted by footy,  <dig> were members of the m <dig> training set,  <dig>  were co-predicted by m <dig> and nnpp <dig> , and  <dig>  were predicted by all four algorithms.
characteristics of the m <dig> genome-wide prediction can be summarized by looking at all  <dig> predictions, or by looking at the  <dig> unique genes, and selecting the predictions closest to the tls.
the two views produce similar results.
approximately 64% of predicted promoters are completely contained in non-coding upstream regions, 50% are on the positive strand, and time of activation distributes as follows: 5% hour  <dig>  23% hour  <dig>  51% hour  <dig>  20% hour  <dig> and 2% hour  <dig>  the strand and hour distributions for all  <dig> genes in the genome are equivalent to the predicted promoter distributions, indicating that there is no strand or temporal preference for the predicted ct σ <dig> promoters.
discussion
the final model produced by the iterative strategy was generated by a training set with three of the original members, ct <dig>  ct681a and ct681b, removed.
an explanation of how these three sequences differ from the remainder would be informative.
the last column of table  <dig> reports that ct <dig> and ct <dig> are both expressed at  <dig> h pi, classifying them as mid-cycle genes.
niehus et al   <cit>  recently demonstrated that chlamydial promoters show a differential response to changes in dna supercoiling that correlates with the lifecycle expression pattern.
specifically, two mid-cycle genes  responded to supercoiling, while three late-cycle genes  did not.
their experimental set included ompa/ct <dig> in the mid-cycle group and omca/ct <dig>  hcta/ct <dig> &ltub/ct <dig> in the late-cycle group.
thus, it is likely that there exists a set of mid-cycle promoters that differ topologically from other promoters to enhance their ability to respond to supercoiling, and this may explain the anomolous characteristics of these promoters that we observed.
a possible explanation for the large number of genes without promoter predictions by the m <dig> model is heterogeneity requiring different models, for example for response to supercoiling.
while investigating the initial model m <dig>  we explored stepwise nominal regression, which allows for the discovery of more than two dependent variable categories.
however, we did not find that a third category was substantiated.
nonetheless, we suspect that future promoter identifications may confirm the existence of more than two promoter populations for σ <dig> in chlamydiales.
a chief limitation of our study includes the challenge of collecting a reliable training set that was discussed earlier.
we also feel that it would be advisable in future studies to relax the range of possible spacer lengths in the duration hmm for increased generalization, which might have allowed the discovery of more promoters in the whole genome analysis.
additionally, it is quite possible that there are structural features downstream from the tls, as well as upstream, which would aid in promoter discovery.
future modeling efforts should extend the region of interest to  <dig> nt downstream from the tls.
the high priority assigned to the duration hmm scores by the sblr procedure reinforces that the duration hidden markov model is an encouraging approach for modelling core promoters, that deserves further development.
also by implementing our model in hmmer our duration hmm is reusable, generalizable, easily adapted to other organisms and open-source.
this approach explicitly incorporates spacing preferences of elements in a likelihood framework.
two natural further developments of this approach would include further iteration of the model development in chlamydia using an expanded training set, exploiting computational criteria and measurements to define expanded training sets.
another possible extension would be to model extended promoter elements using further elaborations of the hidden markov modeling framework.
the ct genome-wide promoter predictions and co-predictions with other algorithms provide the basis for future research in promoter identification.
the fact that 20% of m <dig> predicted promoters were co-predicted by nnpp <dig>  and tss-predict supports the validity of all three predictions.
the expected confirmation of these promoters will augment the list of verified promoters.
however, confirming or rejecting the predictions made by only m <dig> will provide more valuable information.
confirmation will strengthen the current model in a direction that diverges from e. coli, while rejection will add new non-promoter observations that differ from the current training set.
CONCLUSIONS
models m <dig> and m <dig> support the conjecture that measures of dna biophysical criteria along with measures of rnap σ-factor/dna binding collaboratively contribute to a sequence's ability to promote transcription.
whereas a measure of rnap σ-factor/dna binding ensures a sensitive prediction, adding measures of position relative to the tls, stability, curvature, sidd and twist provide specificity.
the stratified k-fold cross-validation of m <dig> indicates that the model performs well by absolute criteria as well as compared to other predictive algorithms.
additionally, there is considerable overlap between the genome-wide predictions of m <dig> and nnpp <dig> , tss-predict and footy.
the modeling procedure we describe here seems especially appropriate for bacterial species where the set of known promoters is limited and the genome is relatively small.
outlook
the model derived by the method described here is a first pass model that serves as proof of concept.
the ct genome-wide promoter predictions, along with co-predictions by nnpp <dig> , tss-predict and footy, will allow researchers to select optimal candidates for validation mapping of transcript 5' ends by primer extension.
as more chlamydial promoters are identified, the model will be updated, and a refined list of promoter predictions may be developed.
more interactions among predictor variables may also be explored.
a final model will provide insight into the process of chlamydial transcription initiation.
then, too, it will be possible to determine if chlamydial promoters differ significantly from those of other bacteria.
authors' contributions
rm participated in conceiving the study, designed the strategies, retrieved data from various websites, conducted the data analysis, performed calculations and wrote the r scripts.
do participated in conceiving the study and provided bacteriological expertise.
da participated in conceiving the study and provided modeling expertise.
research was performed under the advice and supervision of do and da.
all authors contributed to the draft of the paper, and all authors read and approved the final manuscript.
supplementary material
additional file 1
input data file for durahmmer used to build m <dig> duration hmm.
click here for file
additional file 2
parameters of the m <dig> duration hmm that were provided for hmmsearch by durahmmer, an hmmer  <dig> .
<dig> model file.
click here for file
additional file 3
data file used to build m <dig>
click here for file
additional file 4
ct promoters predicted by m <dig>
click here for file
additional file 5
ct promoters predicted by m <dig> and nnpp <dig> .
click here for file
additional file 6
ct promoters predicted by m <dig> and tss-predict.
click here for file
additional file 7
ct promoters predicted by m <dig>  nnpp <dig>  and tss-predict.
click here for file

 identifying mechanisms that regulate gene expression in bacteria is essential for understanding and eventually controlling their pathogenicity. all known bacteria share a well conserved transcriptional holoenzyme, rna polymerase. the rnap is comprised of a 3-subunit catalytic core plus a variable σ-factor subunit that provides dna binding specificity. one of these σ-factors, σ in escherichia coli, participates in the transcription of a majority of genes including those with housekeeping functions. e. coli is the best studied bacterial model with regard to promoter identification and prediction. as such, most promoter predictions are based upon the analysis of e. coli σ promoter data. the earliest collections of e. coli σ promoters revealed the - and - hexamer consensus motifs, ttgaca and tataat, that serve as recognition sites for the and domains of σ. position weight matrices were the first models to quantify the hexamer motifs. pwm models were expanded to quantify the variable-length spacer region between hexamers, which is important in orienting the hexameric motifs for interaction with the sigma binding factors. challenges encountered by pwm models include defining thresholds that are sensitive enough to include known promoters without predicting numerous false positives. most of the quantitative modeling efforts that ensued require training sets comprised of both positive and negative sequences. artificial neural networks have been trained on sequences of identified e. coli promoters and non-promoters. a hidden layer in the ann architecture quantifies interactions among pairs and triplets of nucleotides. the resulting ann scans and scores overlapping sequences, and reports a score in the range that indicates the likelihood of the sequence being a promoter. a time-delay neural network can combine two simple anns with a variable-length spacer region. burden et al measured the distance from the transcription start site to the translation start site of e. coli promoters. they showed that the distribution peaks sharply around nt, and that combining the tss-tls distribution with the nnpp tdnn significantly enhances the specificity of the prediction. in another machine learning approach that has been applied to model promoters, support vector machines were trained on e. coli promoter sequences of length. although the svm approach has the advantage of comprehensively quantifying the primary structure of the upstream region, it does not examine structures of higher order that motivate our approach. a natural extension to pwms that explicitly models an empirical spacing distribution between motifs is given by duration hidden markov models. here "duration" refers to this explicit representation of a spacer length distribution, as opposed to the geometrically distributed lengths that are expected from components of profile hidden markov models. although the variable-length spacer region between hexamers has been incorporated into promoter modeling and predictors before, none of these earlier efforts have integrated an explicit probabilistic representation of the spacer distribution into a reusable predictor as a duration hmm, which is arguably its most natural representation. on the other hand, while duration hmms have been introduced into genome analysis, they have not to our knowledge been applied to modeling transcriptional or translational signals before. bacteria of the genus chlamydia are obligate intracellular parasites that were genetically isolated from other bacteria nearly a billion years ago when they moved into their intracellular environment. in humans, chlamydia infections are responsible for infertility, blindness, arthritis and cardiovascular disease. because chlamydiae have an intracellular life-cycle, standard genetic techniques are often insufficient to study gene regulation. hence, only to promoters have been experimentally verified. however, with a small genome of only about mbp and genes, chlamydia trachomatis makes a good candidate for in silico analysis. surveys of known bacterial promoters suggest that their structures are relatively diverse. additionally, established ct promoters display obvious differences from the established consensus hexamers of e. coli. although σ the ct analog of e. coli σ has dna binding domains homologous to domains and in σ sequence based phylogenetic analysis of bacterial rnap subunits has shown discernable evolutionary distance between the ct and e. coli rnap subunits. therefore, it is plausible that an organism-specific model is appropriate for ct. phylogenetic footprinting takes advantage of relative conservation of motifs among related species. grech et al developed an algorithm that combines e. coli trained pwms and chlamydial phylogenetic footprinting. ct upstream regions are screened with the pwms and the potential promoter hexamers are filtered with an algorithm that accepts only conserved sequences in a consensus of c. trachomatis, c. pneumoniae and c. caviae. although this is a promising approach, because they used an e. coli trained pwm, their results may be strongly influenced by prior expectations that all bacterial promoters are structured as in e. coli. we believe that more development is needed in ab initio approaches for predicting promoters using sequence information directly from the organism under study in combination with biophysical metrics that derive from known models about the biology of transcription in general. this study aims to develop ct promoter models using only known ct promoters. to do so, it considers dna stability and topological features of the upstream region as well as rnap σ-factor/dna binding. as hertz and stormo aptly wrote "... the polymerase needs to bind the dna, open the dna, initiate transcription, and release the promoter for elongation." the tdnns and svms that consider extended promoter sequences are addressing this issue from a sequence perspective. this study utilizes measures that have been developed to quantify stability and other aspects of dna structure. evidence from the profiling of dna curvature, bendability, twist, stability and propensity for stress-induced destabilization in e. coli, b. subtilis, c. trachomatis, plants and vertebrates suggests that there are peaks for these measures near the tss. here we use a stochastic model building procedure that allows for the combination of relevant predictive measures selected from rnap σ-factor/dna binding propensity, as quantified by duration hmms, and structural features of the upstream region, as quantified by biophysical metrics. methods stochastic model building stepwise binary logistic regression, as implemented in spss version statistical software, selects an optimal set of independent variables to classify observations into two populations. logistic regression does not assume a linear relationship between the dependent and independent variables, normal distributions, or homoscedasticity. it does, however, assume independence of observations. we address this requirement in a separate section describing the selection of non-redundant observations. the mathematical model fitted by sblr has the form where i is the number of steps, v through vi are the predictor variables selected, and b through bi are coefficients determined by the analysis. u is the logit for the dependent variable, which means that here, the event is class membership. when p denotes the prob, the equation can be rewritten as selecting a cutoff for p, most commonly, converts p into a classifier. when is the probability threshold, eu = and the classification threshold for u is the effectiveness of a model can be evaluated by its ability to correctly classify the training data. the spss sblr analysis procedure provides many user-defined options. we selected the forward conditional stepwise procedure for all analyses. at each step, a score statistic is calculated for each variable excluded from the model. the score statistic is based on maximum-likelihood estimation criteria and is asymptotically distributed as a χ variable. the variable with the highest significant χ value is entered into the model. if no significant variables remain, then the procedure stops with the current model. similarly there is a mechanism for stepwise removal. after a new model has been generated, score statistics are calculated for all variables in the model. if the p-value for any variable in the model is greater than the probability for stepwise removal, then the variable is removed from the model. we retained the default probabilities for stepwise entry and removal, thus ensuring that the significance of all model variables is less than. potential observations and dependent variable a significant challenge for bioinformaticians is to model data that has been collected by multiple laboratories using different assays, protocols and equipment. this phenomenon is compounded in the study of ct where the organism is metabolically active only inside an infected host-cell. one way to minimize the use of conflicting and/or controversial data is to rely upon reviews written by informed biologists. for this reason, we consulted the reviews of mathews & timms and tan to compile a list of experimentally verified σ promoters. to this list we added promoters that were experimentally verified by grech et al and hefty et al after the previously cited reviews were written. for the purposes of this study, we consider these sequences to be the known ct σ promoters. ant distance from the lead nt of the - hexamer to the tls. breferences: m: mathews & timms; g: grech et al; h: hefty et al; t: tan chour post infection of transcriptional activation files containing the ct genome and genome table were retrieved from the ncbi website,, on july and last modified by ncbi on january an r script was written to extract the nt upstream regions of all protein-coding genes as annotated in the genome table. the regions were verified at several stages throughout the research with other sources, including a ct genome database previously available from los alamos national laboratories and comparable prediction algorithms. table displays the distance from promoter to tls for all training set genes. since the maximum distance is nt, nt was set as the upper limit for data analysis. then, the upstream region was defined as nt to allow for biophysical structures nt upstream from a predicted promoter. for each of the genes listed in table the nt upstream region was parsed into overlapping sliding window sequences of length and step-size each subsequence was labeled according to its parent gene and position occupied in the upstream region: e.g. the first seq was labeled ct046_ because the initial nt is found nts upstream from the ct tls. the dependent variable, promoter, was assigned a if a promoter sequence listed in table was totally contained in seq and otherwise. thus, 1's identify potential promoter observations and 0's identify potential non-promoter observations. cases with upstream positions ≥ and ≤ were selected as potential observations to restrict the analysis to the range of the training set data. the upper bound is nt upstream of the furthest upstream training set promoter and the lower bound is equal to the furthest downstream training set promoter. independent variables the primary variable for promoter prediction is the pattern that characterizes the binding between the rnap σ-factor and dna. here we use duration hmms to describe and quantify rnap-σ66/dna binding. after a set of known promoters is used to train a duration hmm, the duration hmm scans a new sequence to identify the hexamer-spacer-hexamer subsequence that scores the highest with regard to potential binding. the variable hmm_score is assigned the score associated with the highest scoring subsequence, while the variable start denotes the position of the lead nucleotide in the - hexamer and end denotes the position of the last nucleotide of the - hexamer. specifically, a training set of promoter sequences was placed in the file ts.txt. the initial ts.txt contained the contents of table columns "- hex", spacer and "- hex". this file was supplied as input to durahmmer which was used to create a duration hmm with the command: durahmmer - - -s -s -p -u::: -c ts.txt > ts.hmm. the options to the command specify the following model parameters: matched states at the 5' and 3' sequence ends; minimum and maximum spacer lengths of and respectively; a background compositional model of % a, % c, % g, and % t; and spacers should be modeled to have their empirical composition in the training set. the program durahmmer produces a valid hmmer. model file representing a duration hmm. for the final model of this study, the model file and the input data file are provided as additional files and all seq observations from the genes were placed in the file all.txt so that optimal promoters and hmm scores could be calculated by hmmsearch with the command: hmmsearch -e ts.hmm all.txt. we ran hmmsearch with a high e-value because we were interested in combining the score of the maximum scoring hit with other metrics in a composite procedure regardless of its magnitude. in combination with the duration hmm model score described above, we also used the following biophysical metrics of promoter position and structure as possible independent variables for the sblr model: position, which indicates the location of seq in the upstream region relative to the tls. for ct046_ position = measures of curvature and %gc content for each nt upstream region, which were determined by the online bend.it server with a window-size of free energy change of dna melting, bendability and twist angle, which were determined for each nt upstream region by the online plot.it server. all measurements were then averaged over each seq Δg always has a negative sign and is interpreted as greater values having lower stability. for statistical analysis this variable was transformed by stable = -Δg so that the sign is always positive and the interpretation is that larger values have greater stability. stability is also of interest in the immediate downstream region, so positions 27- and 1- were quantified. since the bendability measure increases with rigidity, it was renamed rigid. the twist angle measurement, twist, was not transformed. possible times of expression onset include and hours post infection. mutually exclusive binary variables h h h h h and h were created to mark time of expression onset. stress-induced dna duplex destabilization quantification utilizes structural and energetic properties of dna to measure the propensity for strand separation under negative superhelical stress. a low sidd score indicates a high propensity for strand separation. sidd measurements were determined by the websidd server with default parameters except for open region size = because niehaus et al have shown a time dependent response to chlamydial dna supercoiling, interactions between the time of expression onset and sidd were included. the sidd/hour of onset interaction is quantified by sidd_h# = sidd*h#. for variables based on the entire seq lagged variables were created for the four non-overlapping upstream subsequences of length 32: e.g. for ct046_ curve_l was set equal to the curve value of ct046_ curve_l was set equal to the curve value of ct046_164; curve_l was set equal to the curve value of ct046_196; and curve_l was set equal to the curve value of ct046_ selection of non-redundant observations from potential observations as mentioned earlier, sblr assumes independent observations. to address this requirement, we select for analysis a subset of the overlapping potential observations that are non-redundant with respect to the pair of hexamers that are most likely to bind the rnap σ-factor. if we select only those cases where end = we eliminate all of the redundant optimal hmm hexamer pairs while retaining most optimal hmm instances. table demonstrates how this selection ensures that neighboring optimal hmm instances that match are included only once. six potential observations, ct046_ through ct046_ all contain the verified promoter with hexamer pair tggtta and taaaaa. consequently, they all have promoter = and hmm_score = -. but only ct046_ has end = and is selected to represent the verified ct promoter. similarly, only ct046_ represents the maximal non-promoter hexamer pair ttgtgt and aaaagt with score = -. this process incidentally aligns each selected seq such that the optimal downstream hexamer is at the far right end. this selection process does not eliminate overlapping sequences, but it does eliminate overlapping likely binding sites. ct046_ and ct046_ overlap a great deal. however, the last hexamer of ct_046_ is not present in ct046_ and the first hexamer of ct_046_ does not appear in ct046_ it should be noted that although each training set gene begins with the same number of potential observations, this selection process causes the number of selected non-redundant observations to differ among genes. each gene starts with around potential observations with promoter = for each verified promoter, and around 325-40- = potential observations with promoter = however, selection for non-redundant observations always results in the number of designated non-promoters being reduced to approximately while selecting sequences with non-redundant hmm_scores does mitigate the problem of dependent observations, it may not entirely eliminate it. while there are numerous studies that affirm the robustness of baysian discriminant analysis with regard to violating the assumptions of a linear relationship between the dependent and independent variables, normal distributions, and homoscedasticity, we could not find similar studies regarding the robustness of logistic regression. an alternative to the current analysis would be to use stepwise discriminant analysis, knowing that we are violating some assumptions. there are versions of logistic regression, including generalized estimating equations, that are specifically designed for correlated data such as longitudinal studies. in these procedures there are subject variables and within-subject variables. it might be possible to force this study data into such a format, but as yet there are no readily available stepwise procedures to scan multiple possible predictors. a final alternative would be to select non-overlapping sequences with the penalty of losing information and perhaps introducing a selection bias. sblr is a procedure for model identification. it is only after a model has been identified that it can be evaluated for independence. given that, we elected to analyze the non-redundant observations with sblr and then examine the error terms for independence. in time series analysis, this is done by checking that the error term is normally distributed with zero mean, and that autocorrelations and partial autocorrelations of the error term are not significant. iterative modeling strategy sources of error that could lead to misclassification include imprecise laboratory procedures in defining and identifying promoters, presence of more than one promoter population, failure to include relevant predictor variables, and random variation. to minimize the first two error sources, an iterative strategy was developed. duration hmm iteration addresses error source, while sblr iteration addresses source. duration hmm iteration minor modifications in the configuration of the training set promoters can improve classification accuracy. to accomplish this, we allowed each promoter to vary within a neighborhood that extends the sequence by nts on each side. a limit of nts ensures that a modified hexamer will not locate completely outside of the original promoter sequence. for example, when the promoter ct is extended, it becomes ttgtttgcagagtttttattttaaatatgttataatctgtc, with the bolded nts marking the extensions. initially, a duration hmm is determined by the original, non-extended, promoter training set. then the set of extended promoters is searched for the highest scoring instance of the duration hmm in each extended sequence. if a high-scorer is not the same as the original promoter, it replaces the original in the training set. the iteration continues until stabilization. for the final model, ct was modified to ttgcagagtttttattttaaatatgttataat. sblr iteration deletion and subsequent replacement of members of the training set can eliminate promoters that are likely to be members of a different promoter population. this is accomplished via the iterative scheme diagrammed in figure initially, the complete set of verified promoters determines the duration hmm and the independent observations selected for sblr analysis. sblr delivers a mathematical model that produces a predicted probability of class membership for each observation. a threshold on p of. is used to classify each observation as a predicted promoter or non-promoter. for those cases where promoter = we also use the value of p to determine when a promoter appears to be an outlier and should be eliminated from the training set. after observing the probabilities, a retention threshold on p between and. is established. if a training gene has only one identified promoter and that promoter has a p less than the retention threshold, then all observations for that gene are deleted from the analysis. similarly, if a training set gene has two identified promoters and they are both selected for deletion, all observations for that gene are deleted. however, if a training set gene has two identified promoters and only one is selected for deletion, all upstream observations for that gene remain in the analysis dataset and only observations within the remaining promoter are assigned promoter = modifying the training set in any way necessitates the determination of a new duration hmm, which in turn determines which observations will be aligned such that end = and subsequently included in the next sblr analysis. the iteration process continues until the training set stabilizes. stratified k-fold cross-validation once the final training set and model are selected, it is necessary to validate the model to ensure against over-fitting and to allow for comparisons with algorithms trained on other datasets. in the case of dichotomous classification, stratified k-fold cross-validation partitions the training set into k subsamples such that each subsample has approximately the same proportions of class membership. here we designate each training gene as a subsample; hence k equals the number of genes in the training set. then, one gene is retained as a validation set while the remaining genes are used as training data. evaluation measures are calculated by aggregating the results of each validation set. comparable algorithms the following three algorithms were used to compare performance and to identify co-predictions with the model developed in this study: nnpp, tss-predict, and footy. nnpp is an online time-delay neural network that is accessible for promoter predictions at. we used the following options: organism = prokaryote and minimum promoter score = to define promoters in the nt upstream region of all ct genes. for the support vector machine algorithm tss-predict, the top two ranking predictions for each ct gene are posted as supplementary material at doi: /j.combiolchem.. the ct promoters predicted by footy, an algorithm that utilizes phylogenetic footprinting, are reported directly in the publication that describes the algorithm. r scripts scanned the promoters predicted by nnpp and tss-predict for matches with the promoters predicted by the study model. an nnpp match was declared when the study prediction was contained within the nt nnpp prediction. a tss_predict match was declared when the tss_predict predicted hexamer pair was contained within the study prediction.  finding the best model the initial model, m utilizes the initial training set of promoters with observations from their parent genes. the duration hmm model converged after one iteration, modifying the alignment of promoters. for all models, table reports the variables that were selected for the model and evaluation measures. if tp = true positive, fp = false positive, tn = true negative and fn = false negative, then sensitivity or recall = tp/, specificity = tn/, positive predictive value or precision = tp/, negative predictive value = tn/, and accuracy = /. the total number of observations for each model differs according to the promoter training set being used. training set ct681a ct681b ct681a ct681a +stable1_37 -position +curve_l32 -gc_l128 +rigid_l96 +stable1_37 -gc_l32 -position +curve_l32 -curve_l64 -gc_l128 +stable1_37 -position +curve_l32 -stable_l32 +sidd_h24 -curve_l128 -sidd_l128 +stable1_37 -position +curve_l32 -stable_l32 -stable27_37 ppv or athe variables are listed in order of entrance into the model and the sign indicates the sign of the coefficient. broc analysis area under the curve for model m of the promoters were classified correctly, with false positives. there is always the possibility that these are yet to be recognized promoters, but at this point they are counted as misclassifications. for the verified promoters that were missed, the predicted probabilities ranged from to. since a natural separation appeared to between and, p = was selected as the retention threshold and promoters ct ct681a, ct681b and ct were deleted from the training set for the next model, m the duration hmm model for m converged after two iterations, modifying the alignment of promoters. table shows that m classified the modified training set perfectly, indicating that perhaps too many promoters had been deleted from the original training set. the retention threshold was reset to and ct was reinstated for model m the duration hmm model for m converged after one iteration. table displays the alignments of the promoters that were modified. m also classified the modified training set perfectly. again the results indicated that the next model, m should reset the retention threshold to and reinstate ct however, table reports that m is not as good as models m and m because of classification errors. given two models, one training set a subset of the other, that both classify their respective training sets with 100% accuracy, we reasoned that the model trained on the largest set would provide the most sensitive genome-wide prediction. thus, m was selected as the best and final model because of the perfect classification with the largest training set. the complete data file used to build m additional file is supplied so that others may replicate or modify the model. finally, the error terms of m were checked for independence. residuals, promoter - p, were calculated for all selected observations and shown to be normally distributed with zero mean. additionally, the autocorrelations and partial autocorrelations of the residuals were not significant. thus, the independence assumption of sblr was not violated by this model. aggregated results of the stratified k-fold m cross-validation are reported in the last column of table for the genes in the m training set, promoters were not identified and there were false-positive predictions. the incorrect classifications are most likely due to incomplete representation of the sample space, but may indicate additional populations or absent predictors. model interpretation the m duration hmm describes and quantifies the rnap-σ66/dna binding observed in the training set. a visualization of the m parameters is shown in figure the - hexamer is dominated by the initial ttg motif, while the initial t with frequent as and ts describe the - hexamer. the c and g compositions of the spacer region are much smaller than those of the genome. spacer lengths of predominate, while spacers of length are absent. the input data file for durahmmer and the resulting output data file are provided as additional files and the output data file is an hmmer. model file which supplies the parameters of the m duration hmm to hmmsearch. complete documentation for the contents of the file can be found in the hmmer user's guide at. briefly, the first lines are header information with the main model section following. there are lines for each of the possible nodes. the first and last nodes refer to the - and - hexamers, while nodes through refer to possible spacer positions. the first line for each node displays the contribution to the final score for the corresponding nucleotide matching a, c, g or t. the third line is particularly relevant to nodes through which correspond to spacer nucleotides through as nucleotides in these positions may or may not be present in the sequence being scored due to variable spacer length, the third line provides the odds of transitioning to another spacer nucleotide or to the - hexamer. the m prediction equation generated by sblr is: being the strongest predictor, hmm_score is selected in the first step of the sblr procedure. the prediction equation for step one is using a classification cutoff of p = and setting u = yields hmm_score = as the threshold for step classification. at step 14/ promoters and 2220/ non-promoters were classified correctly. thus, the remaining eight model variables moved promoters with hmm_score <. to promoter classification and non-promoters with hmm_score ≥. to non-promoter classification. the predictor variables and their coefficients describe the verified promoters and their upstream regions. promoters have high hmm_score and low position. the near upstream region is curved and unstable, whereas the further upstream region is uncurved and unstable under superhelical stress. for late-cycle genes where expression onset occurs at h pi, the effect of superhelical stress is less than at other times. the upstream characteristics may reflect transcription factor binding and/or additional interaction with the rnap holoenzyme. the interpretation of the positive coefficient for stable1_ is more subtle. in the second step of the sblr, four observations change from fp to tn and observations change from fn to tp. the means of stable, stable1_ and stable33_ are all larger in the second group than in the first. although stable33_ shows the greatest mean difference, the most statistically significant is stable1_ model exercise: predicting promoters for the ct genome finally, the m model was used to predict promoters for the entire ct genome. additional file reports predicted promoters in unique genes, along with their hmm scores and genome locations. thus, for of the total ct genes, this model does not find any 32-mers with a probability >. this suggests a conservative prediction that emphasizes specificity over sensitivity. other explanatory factors may include alternate binding patterns for σ alternative σ-factors, and operon configurations. there was a substantial overlap among predictions by different methods. additional file lists the promoters co-predicted by m and nnpp, while additional file lists the promoters co-predicted by m and tss-predict. additional file reports the promoters co-predicted by m nnpp and tss-predict. all predictions are for = position = consistent with the range of the modeling procedure. of the promoters predicted by footy, were members of the m training set, were co-predicted by m and nnpp, and were predicted by all four algorithms. characteristics of the m genome-wide prediction can be summarized by looking at all predictions, or by looking at the unique genes, and selecting the predictions closest to the tls. the two views produce similar results. approximately 64% of predicted promoters are completely contained in non-coding upstream regions, 50% are on the positive strand, and time of activation distributes as follows: 5% hour 23% hour 51% hour 20% hour and 2% hour the strand and hour distributions for all genes in the genome are equivalent to the predicted promoter distributions, indicating that there is no strand or temporal preference for the predicted ct σ promoters. discussion the final model produced by the iterative strategy was generated by a training set with three of the original members, ct ct681a and ct681b, removed. an explanation of how these three sequences differ from the remainder would be informative. the last column of table reports that ct and ct are both expressed at h pi, classifying them as mid-cycle genes. niehus et al recently demonstrated that chlamydial promoters show a differential response to changes in dna supercoiling that correlates with the lifecycle expression pattern. specifically, two mid-cycle genes responded to supercoiling, while three late-cycle genes did not. their experimental set included ompa/ct in the mid-cycle group and omca/ct hcta/ct &ltub/ct in the late-cycle group. thus, it is likely that there exists a set of mid-cycle promoters that differ topologically from other promoters to enhance their ability to respond to supercoiling, and this may explain the anomolous characteristics of these promoters that we observed. a possible explanation for the large number of genes without promoter predictions by the m model is heterogeneity requiring different models, for example for response to supercoiling. while investigating the initial model m we explored stepwise nominal regression, which allows for the discovery of more than two dependent variable categories. however, we did not find that a third category was substantiated. nonetheless, we suspect that future promoter identifications may confirm the existence of more than two promoter populations for σ in chlamydiales. a chief limitation of our study includes the challenge of collecting a reliable training set that was discussed earlier. we also feel that it would be advisable in future studies to relax the range of possible spacer lengths in the duration hmm for increased generalization, which might have allowed the discovery of more promoters in the whole genome analysis. additionally, it is quite possible that there are structural features downstream from the tls, as well as upstream, which would aid in promoter discovery. future modeling efforts should extend the region of interest to nt downstream from the tls. the high priority assigned to the duration hmm scores by the sblr procedure reinforces that the duration hidden markov model is an encouraging approach for modelling core promoters, that deserves further development. also by implementing our model in hmmer our duration hmm is reusable, generalizable, easily adapted to other organisms and open-source. this approach explicitly incorporates spacing preferences of elements in a likelihood framework. two natural further developments of this approach would include further iteration of the model development in chlamydia using an expanded training set, exploiting computational criteria and measurements to define expanded training sets. another possible extension would be to model extended promoter elements using further elaborations of the hidden markov modeling framework. the ct genome-wide promoter predictions and co-predictions with other algorithms provide the basis for future research in promoter identification. the fact that 20% of m predicted promoters were co-predicted by nnpp and tss-predict supports the validity of all three predictions. the expected confirmation of these promoters will augment the list of verified promoters. however, confirming or rejecting the predictions made by only m will provide more valuable information. confirmation will strengthen the current model in a direction that diverges from e. coli, while rejection will add new non-promoter observations that differ from the current training set.  models m and m support the conjecture that measures of dna biophysical criteria along with measures of rnap σ-factor/dna binding collaboratively contribute to a sequence's ability to promote transcription. whereas a measure of rnap σ-factor/dna binding ensures a sensitive prediction, adding measures of position relative to the tls, stability, curvature, sidd and twist provide specificity. the stratified k-fold cross-validation of m indicates that the model performs well by absolute criteria as well as compared to other predictive algorithms. additionally, there is considerable overlap between the genome-wide predictions of m and nnpp, tss-predict and footy. the modeling procedure we describe here seems especially appropriate for bacterial species where the set of known promoters is limited and the genome is relatively small. outlook the model derived by the method described here is a first pass model that serves as proof of concept. the ct genome-wide promoter predictions, along with co-predictions by nnpp, tss-predict and footy, will allow researchers to select optimal candidates for validation mapping of transcript 5' ends by primer extension. as more chlamydial promoters are identified, the model will be updated, and a refined list of promoter predictions may be developed. more interactions among predictor variables may also be explored. a final model will provide insight into the process of chlamydial transcription initiation. then, too, it will be possible to determine if chlamydial promoters differ significantly from those of other bacteria. authors' contributions rm participated in conceiving the study, designed the strategies, retrieved data from various websites, conducted the data analysis, performed calculations and wrote the r scripts. do participated in conceiving the study and provided bacteriological expertise. da participated in conceiving the study and provided modeling expertise. research was performed under the advice and supervision of do and da. all authors contributed to the draft of the paper, and all authors read and approved the final manuscript. supplementary material additional file 1 input data file for durahmmer used to build m duration hmm. click here for file additional file 2 parameters of the m duration hmm that were provided for hmmsearch by durahmmer, an hmmer. model file. click here for file additional file 3 data file used to build m click here for file additional file 4 ct promoters predicted by m click here for file additional file 5 ct promoters predicted by m and nnpp. click here for file additional file 6 ct promoters predicted by m and tss-predict. click here for file additional file 7 ct promoters predicted by m nnpp and tss-predict. click here for file
model based genome-wide predictions are provided so that optimal promoter candidates can be experimentally evaluated, and refined models developed. this work provides a baseline model that can evolve as new chlamydia trachomatis σ promoters are identified with assistance from the provided genome-wide predictions. although it has been shown that transcriptional mechanisms are similar among various bacteria, it is quite possible that the differences between escherichia coli and chlamydia trachomatis are large enough to recommend an organism-specific modeling effort. here we present an iterative stochastic model building procedure that combines such biophysical metrics as dna stability, curvature, twist and stress-induced dna duplex destabilization along with duration hidden markov model parameters to model chlamydia trachomatis σ promoters from experimentally verified sequences. initially, iterative duration hidden markov modeling of the training set sequences provides a scoring algorithm for chlamydia trachomatis rna polymerase σ66/dna binding. in addition, most promoter models have been trained on data from escherichia coli. the resulting model predicts the final training set with a high degree of accuracy and provides insights into the structure of the promoter region. promoter identification is a first step in the quest to explain gene regulation in bacteria.
 promoter identification is a first step in the quest to explain gene regulation in bacteria. it has been demonstrated that the initiation of bacterial transcription depends upon the stability and topology of dna in the promoter region as well as the binding affinity between the rna polymerase σ-factor and promoter. however, promoter prediction algorithms to date have not explicitly used an ensemble of these factors as predictors. in addition, most promoter models have been trained on data from escherichia coli. although it has been shown that transcriptional mechanisms are similar among various bacteria, it is quite possible that the differences between escherichia coli and chlamydia trachomatis are large enough to recommend an organism-specific modeling effort.  here we present an iterative stochastic model building procedure that combines such biophysical metrics as dna stability, curvature, twist and stress-induced dna duplex destabilization along with duration hidden markov model parameters to model chlamydia trachomatis σ promoters from experimentally verified sequences. initially, iterative duration hidden markov modeling of the training set sequences provides a scoring algorithm for chlamydia trachomatis rna polymerase σ66/dna binding. subsequently, an iterative application of stepwise binary logistic regression selects multiple promoter predictors and deletes/replaces training set sequences to determine an optimal training set. the resulting model predicts the final training set with a high degree of accuracy and provides insights into the structure of the promoter region. model based genome-wide predictions are provided so that optimal promoter candidates can be experimentally evaluated, and refined models developed. co-predictions with three other algorithms are also supplied to enhance reliability.  this strategy and resulting model support the conjecture that dna biophysical properties, along with rna polymerase σ-factor/dna binding collaboratively, contribute to a sequence's ability to promote transcription. this work provides a baseline model that can evolve as new chlamydia trachomatis σ promoters are identified with assistance from the provided genome-wide predictions. the proposed methodology is ideal for organisms with few identified promoters and relatively small genomes. 