BACKGROUND
most eukaryotic genomic dna is wrapped in nucleosomes, which occlude and strongly distort the wrapped dna.
accumulating evidence shows that the dna sequence itself is highly predictive of nucleosome positioning in vivo  <cit> , and that nucleosome positioning is closely related to chromosome functions  <cit> .
a fast software tool for predicting nucleosome positioning is highly desirable.
several statistical methods for nucleosome positioning prediction have been proposed in the literature.
in  <cit>  a method was proposed based on cross-correlation with a nucleosome dna sequence signature.
in  <cit>  a markov model was used together with consideration of steric exclusion and thermodynamic equilibrium.
in  <cit> , a support vector machine  was trained based on the differential k-mer usage between nucleosome and linker dnas.
in  <cit> , the authors proposed an n-score model to discriminate nucleosome and linker dnas using wavelet energies as covariates in a logistic regression model.
in  <cit> , a web-interface called "nuscore" for estimation of the affinity of histone core to dna and prediction of nucleosome positioning was developed based on the dna deformation energy score.
in  <cit> , the model from  <cit>  was improved by incorporation of differential k-mer usage  tracts, which are strongly disfavored by nucleosomes).
this model can be further improved by accounting for nucleosome-nucleosome interaction  <cit> .
while the nucleosomal features are universal, eukaryotic genomes vary in nucleosomal repeat length  <cit>  and base composition.
the nucleosomal repeat length is dictated by the length distribution of linker dnas that separate neighboring nucleosomes, and it determines the overall nucleosome density in the chromatin fiber.
the contribution of this paper is a duration hidden markov model and a software tool called nupop for genome-wide nucleosome positioning prediction.
we show that incorporation of linker length information can achieve better sensitivity in prediction.
in addition, we propose a re-scaling method to adjust for base composition variation when using yeast models to make predictions for other species.
a relatively superior performance of this approach is established by comparing it with other existing tools.
methods
model
the hidden markov model  has been known for decades.
an excellent and famous tutorial is rabiner's  <dig> paper  <cit> , in which the model, algorithms, and applications were carefully and thoroughly reviewed.
a conventional hmm implicitly assumes a geometric duration distribution for each state, which can be wrong in real applications.
modeling the explicit duration of each state can improve the prediction of hmms .
we model each chromosomal dna sequence with a duration hidden markov model  of two oscillating states: nucleosome  and linker , where the nucleosome state has a fixed length of  <dig> bp, and the linker state has a variable length.
we assume that at the end of each state, the chain must transit to the other state; additionally, a complete chromatin sequence must start with and end in a linker state.
we trained a 4th order time-dependent markov chain for the the n state, and a homogeneous 4th order markov chain for the l state to distinguish the k-mer usage preferences for k up to  <dig> between the nucleosome and linker states as shown in other methods, e.g.,  <cit>  .
let e = e <dig>  ..., e <dig> be a nucleosome dna sequence.
let pn be the probability of observing e as a nucleosome, computed as the product of probabilities for both watson and crick strands under the 4th order markov chain model.
we assume that the linker dna length of a given species has an unknown distribution fl defined for k =  <dig>  ..., τl .
an observed linker dna sequence e = e <dig>  ..., ek carries two pieces of information, the length is k bp, and given which, the emitted letters are e <dig>  ..., ek.
let gl denote the homogeneous markov chain model for the linker dna .
then observing e as a linker dna has probability
suppose x = x <dig>  ..., xn is a genomic dna sequence of length n, where xi = a/c/g/t.
let z = z <dig>  ..., zn be the corresponding hidden state path, where zi =  <dig> if xi is covered by a nucleosome state, and  <dig> otherwise.
suppose that the path z = z <dig>  ..., zn partitions x into k consecutive nucleosome or linker state blocks, in which the nucleosome blocks have a uniform length of  <dig> bp, whereas the length of linker blocks may vary.
we denote these blocks as y = y+, ..., yb, and their state identification as s = s <dig>  ..., sb, where si =  <dig> if yi is nucleosome state, and  <dig> otherwise.
the probability of observing  is given by
where π <dig> and πe stand for the probabilities that the chain initializes and ends with state s <dig> and sk respectively, and i is an indicator function.
since we assume that a complete chromatin sequence must start with and end in a linker state, π <dig> = πe =  <dig>  we define the nucleosome occupancy at a specific position i, denoted oi, as the posterior probability that zi =  <dig>  i.e.,
we also define the histone binding affinity score at position i as the log likelihood ratio for the region xi- <dig>  ... xi, ..., xi+ <dig> to be a nucleosome vs. a linker, i.e.,
given the models pn, gl and fl, the optimal path z can be found by the standard viterbi algorithm, and the nucleosome occupancy score can be estimated using forward and backward algorithms.
data and model training
we utilized the  <dig>  yeast nucleosome dna reads from  <dig> pyrosequencing published in  <cit>  for model training and assessment.
among  <dig>  reads that each were mapped to a unique region of the yeast genome, we first selected reads of length between  <dig> and  <dig> bp.
if multiple such reads existed for the same nucleosome, we selected the one with the highest blast score.
the resulting non-redundant set of  <dig>  nucleosome sequences were center aligned to train the nucleosome model pn.
the 4th order time dependent markov chain can be defined by the base composition at the first position qn, and the transitional probabilities qn, qn, qn, qn, for k =  <dig>  ...,  <dig>  xi = a/c/g/t, i =  <dig>  ...,  <dig>  where the subscript k, i index the positions within a nucleosome.
these probabilities are trained using the corresponding observed fractions or conditional fractions based on the center alignment, with a three bp moving average .
we further identified  <dig>  reads-free regions of length 7- <dig> bp to train the linker state model gl.
the 4th order homogeneous markov model for the linker dnas can be completely defined by the stationary base composition ql, and the transition probabilities ql, ql, ql, ql.
by "homogeneous", we mean that these probabilities are all constants as functions of i. these probabilities were trained using their observed values as in the nucleosome model.
for example, ql was trained by calculating the fraction of occurrences of transitions from any four letters to the fifth letter in the selected putative linker dna sequences.
our initial nucleosome/linker model was trained using the yeast data.
a complication arises when predicting nucleosomes for other species because organisms may differ significantly in their dna base composition.
we propose to scale up or down the probabilities in the markov models by a factor determined by the difference of the base composition between the current species and yeast.
for example, in c. elegans, the fraction of a plus t bases is  <dig>  compared to  <dig>  in yeast.
for a specific transition probability qn at any specific nucleosomal position defined for yeast, we scaled it up as qn ×  <dig> / <dig>  for the corresponding transition probability at that given position for c. elegans.
likewise the transition probabilities for g/c will be scaled down by a factor of  <dig> / <dig> .
all the re-scaled probabilities are then normalized.
the same re-scaling applies to the linker model.
we shall use simulations below to show that re-scaling improves prediction regarding sensitivity and false discovery rate.
using the trained nucleosome model  and linker model , we further train the linker dna length distribution as follows.
we assume that the linker dnas in any given species or cell type have a maximum length τl =  <dig> bp.
this algorithm contains the following steps:
<dig>  initialize the algorithm with a uniform distribution for fl for k =  <dig>  ... τl where τl is the maximum allowable linker length.
<dig>  use the forward and backward algorithm to obtain the posterior expectation of fl for each k. for a sequence x = x <dig>  ..., xn, let nk be the number of linker dnas of length k. then
for k =  <dig>  ..., τl .
for k =  <dig>  ..., τl .
<dig>  update the empirical linker length distribution from step  <dig> using a kernel smoothing method as follows:
where k is the standard gaussian kernel and h is the bandwidth parameter optimally chosen by the leave-one cross-validation method as in  <cit> .
<dig>  use the updated linker length distribution from step  <dig> to compute the nucleosome occupancy and optimal positioning.
compared to viterbi training , using the posterior expectation obtained in eq.
combined with the kernel method in eq.
performs overwhelmingly better in minimizing the summed square errors  .
in the developed software tool nupop, we have trained the linker dna length distributions for  <dig> different species including human, mouse, rat, zebrafish, d. melanogaster, c. elegans, s. cerevisiae, c. albicans, s. pombe; a. thaliana and maize.
the linker dna length distribution  for each species has been trained by scanning the corresponding genome sequences based on τ l =  <dig>  we found that the re-scaled nucleosome and linker profiles, together with the trained linker length distribution, not only roughly recover the genome-wide base compositions, but also the dinucleotide frequencies for different species.
the frequency of each single or di-nucleotide in simulated genomes typically differs by ≤ 1% from that observed in the corresponding real genomes .
as different cell types from the same organism  can exhibit quite different linker dna length distributions  <cit> , a useful future refinement would utilize high quality nucleosome maps for the given cell type, when such data become available.
software tools
we have developed a software tool called nupop, implemented in three different formats: an r package tested for windows xp, linux and mac os x; a stand-alone fortran program; and an nupop web server, all available from http://nucleosome.stats.northwestern.edu.
the r package is built upon the fortran program.
it provides additional handy functions to visualize the resulting viterbi  and nucleosome occupancy predictions.
both the r package and fortran program can handle a genomic sequence of any length with a ram demand < <dig> m bytes.
the predicted results are stored locally in the working directory.
the web server provides an interface through which the user can submit their own sequence up to  <dig> k bp in length for fast online prediction.
when making a prediction, the user is required to specify which species the genomic sequence is from.
if the species is not on the list, nupop will calculate the base composition of the input dna sequence and then choose the nucleosome and linker models from a species that has the most similar base composition.
an alternative model with a 1st order time-dependent markov chain for the nucleosome state and a homogeneous 1st order markov model for the linker state, trained in the same way, is also implemented in nupop as an option.
RESULTS
updating fl improves prediction
updating the linker length distribution not only helps recover the true nucleosome density, but also improves prediction.
we demonstrated this by simulation as follows.
we simulated  <dig> genomic sequences with the 4th and 1st order yeast models respectively, each containing  <dig>  nucleosomes and  <dig>  linkers.
the linker dna length was simulated from a normal distribution  and a gamma distribution .
if a nucleosome is predicted within ± <dig> bp of a true nucleosome, we define it as a correct positive prediction.
the rate of the correct positive prediction, referred to as sensitivity, is defined as the percentage of the  <dig>  nucleosomes that are correctly predicted.
in addition, we include the false discovery rate , defined as the fraction of the predicted nucleosomes that reside > ± <dig> bp away from any true nucleosomes, as the second measure for model performance.
in analogy to statistical hypothesis testing, the sensitivity measures the power of prediction, while the fdr measures the fraction of type i errors in the positive claims.
the results are presented in table  <dig> and  <dig>  in both cases, the linker length was initialized as a uniform distribution with τl =  <dig>  compared to the gamma distribution, the normal distribution is relatively flatter.
thereby updating the linker length distribution did not significantly change the total number of predicted nucleosomes .
the sensitivity increased on average by ~ 4-5% and the fdr dropped by ~ 5% after one update .
further updating continued to improve the prediction until it stabilized after four iterations.
in contrast, the gamma model is much more skewed.
the uniform linker length distribution resulted in an under-estimation of the total number of nucleosomes.
by four updatings, the sensitivity increased by 8%, while fdr remained relatively more stable.
total predictions, sensitivity, and false discovery rate  are the averages  based on  <dig> repeated simulations.
for each simulation a genomic sequence consisting of  <dig> nucleosomes and  <dig> linkers were simulated using the 1st and 4th order yeast models.
the linker length distribution was initialized as uniform on  <dig> ...,  <dig>  and was iteratively updated in the dhmm.
results are shown after  <dig>   <dig>  and  <dig> updates.
total predictions, sensitivity, and false discovery rate  are the averages  based on  <dig> repeated simulations.
for each simulation a genomic sequence consisting of  <dig> nucleosomes and  <dig> linkers were simulated using the 1st and 4th order yeast models.
the linker length distribution was initialized as uniform on  <dig> ...,  <dig>  and was iteratively updated in the dhmm.
results are shown after  <dig>   <dig>  and  <dig> updates.
we also observe that under the same setting and condition, the fourth order model performs slightly but uniformly better than the first order model in both sensitivity and fdr.
this is given that the true models are known and we scan the sequence using the true models.
in theory, the first order markov chain model is nested in the fourth order model.
therefore if the true model is the first order, a well trained fourth order model will have the same prediction power as the first order model, but not vice versa.
since training a higher order markov chain model requires more data, inadequate training can undermine the prediction power.
re-scaling vs. not re-scaling
to illustrate the advantages of re-scaling, we re-scaled the yeast profiles according to the base composition of the maize genome .
using the re-scaled profiles we simulated  <dig> genomic sequences that each contain  <dig>  nucleosomes and  <dig>  linkers.
the linker dna length followed the same two distributions as in table  <dig> and  <dig>  we compare the prediction results from the scaled and non-scaled yeast profiles in table  <dig> and  <dig>  we found using the re-scaled  profile yields a lower fdr than using the yeast profile.
in addition, updating the linker length under the correct profile consistently improves the sensitivity and fdr until prediction stabilizes.
in contrast, while using the yeast profile to scan the simulated maize-like genome, the prediction drastically deteriorates as the linker length updating proceeds.
the same simulation was repeated on other species including human and c. elegans, where the base composition is similar to yeast .
unsurprisingly, the results from the scaled yeast profile were still better than those from the original yeast profile in terms of both sensitivity and fdr, while the difference is much smaller than for the maize case .
total predictions, sensitivity, and false discovery rate  are the averages  based on  <dig> repeated simulations.
for each simulation a maize-like genomic sequence consisting of  <dig> nucleosomes and  <dig> linkers were simulated using the re-scaled 1st and 4th order yeast models.
each sequence was scanned using the true models  and the yeast models with an initial uniform linker length distribution on  <dig> ...,  <dig>  the results after  <dig>   <dig>   <dig> updates of linker length distribution are compared.
total predictions, sensitivity, and false discovery rate  are the averages  based on  <dig> repeated simulations.
for each simulation a maize-like genomic sequence consisting of  <dig> nucleosomes and  <dig> linkers were simulated using the re-scaled 1st and 4th order yeast models.
each sequence was scanned using the true models  and the yeast models with an initial uniform linker length distribution on  <dig> ...,  <dig>  the results after  <dig>   <dig>   <dig> updates of linker length distribution are compared.
nupop vs. other software tools
we briefly assess the prediction performance of nupop by comparing it with two existing methods: the n-score method of  <cit>   and the markov model/thermodynamic equilibrium method of  <cit>  .
as the exact genome-wide nucleosome positioning map is unknown, we utilize the  <dig>   <dig> high-throughput sequence reads to identify well-defined nucleosomes.
we first selected sequences of length between 130- <dig> bp and constructed a reads-based occupancy map.
the reads-occupancy score at a specific position is defined as the number of reads that covered this position.
then we calculated the moving average of this reads occupancy score using a  <dig> bp window.
a sharp peak in the average occupancy curve indicates a nucleosome with well-defined positioning.
considering that the average linker dna length is  <dig> bp in yeast  <cit> , we quantified the sharpness of the peak by calculating the slope from the peak point to the up-/down-stream  <dig> bp point on the average occupancy curve.
we set one condition for the peak to be selected as the center of a nucleosome to be that the absolute value of slope from either side should be >  <dig> .
secondly, we required that the peak height itself must be ≥  <dig> , i.e., a well-defined nucleosome must be testified by at least two well-overlapped reads.
we chose the threshold value as  <dig>  instead of  <dig>  because the overlap of the two reads can be less than  <dig> bp, resulting in a peak on the moving average curve slightly lower than  <dig> .
with these criteria, a total of  <dig>  well-defined nucleosomes are selected from the  <dig> chromosomes of yeast.
a snapshot of a region with many selected well-defined nucleosomes is presented in figure  <dig>  figure  <dig> provides a snapshot of nucleosome occupancy predicted by nupop together with the reads-occupancy.
we first assess the sensitivity of predictions from nupop using the well-defined control set.
if there is a predicted nucleosome within ±k bp of any well-defined nucleosomes , we count this as one correct prediction.
we varied k from  <dig>   <dig>  ...  <dig>   <dig> to investigate the sensitivity behavior at different precision thresholds.
the n-score model predicted  <dig>  nucleosomes.
the current software tool for mm/te method does not provide viterbi predictions, but only the nucleosome occupancy scores.
therefore we calculated the moving average of the occupancy score using a  <dig> bp moving window.
the resulting peaks were treated as the centers of predicted nucleosomes.
if two peaks reside within  <dig> bp, we discarded the one with smaller moving average of occupancy score.
this procedure identified  <dig>  predicted nucleosomes .
likewise, using the occupancy scores from nupop, we identified  <dig>  and  <dig>  predicted nucleosomes under the 4th and 1st order models respectively.
in figure 3a, we compare the sensitivity estimates from the 4th order model of nupop with the n-score and mm/te methods at different threshold values of prediction accuracy .
as the sensitivity tends to increase with an increase in the total predictions, we further selected  <dig>  and  <dig>  best predictions from nupop model  to compare with the n-score method and in figure 3b and 3c respectively.
the sensitivity results suggest that the predictions from nupop outperforms the other two methods in two senses.
firstly, the sensitivity from nupop is  <dig> - <dig> % higher than the other two methods at different threshold values .
secondly, while controlling the total predictions to be the same, nupop has ~  <dig> - <dig> % higher sensitivity than the mm/te method , when the precision threshold is ≤ ± <dig>  as the precision threshold gets less stringent, the difference attenuates and eventually vanishes.
the contrast between nupop and the n-score method is even larger as shown in figure 3c.
as a further comparison, we computed the predictions of the dhmm method under a uniform linker length distribution defined on  <dig>   <dig> ...,  <dig>  this method predicted  <dig>  nucleosomes under the 4th order models, achieving a sensitivity  <dig> - <dig> % higher than the mm/te method .
when controlling the total predictions to be the same as mm/te method or n-score method, the resulting sensitivity curve almost perfectly overlaps with that from nupop.
therefore we omitted these results from figure 3b and 3c.
one could further attempt to evaluate the false positive rate , measuring the fraction of linker regions that were falsely classified as nucleosome regions .
this task requires well-defined linker regions.
a problem, however, is that the average length of linker dnas in yeast  is smaller than the dispersion in lengths of the nucleosome dnas as isolated biochemically .
thus existing nucleosome maps lack the precision needed to define such short linker dnas.
moreover, various sampling biases such as the dna sequence preferences of the micrococcal nuclease used to liberate nucleosomes biochemically  could yield longer genomic regions that are free of recovered nucleosome dna reads even if they are actually nucleosome occupied  <cit> .
attempts to evaluate the fpr given these problems in the data could result in misleading conclusions.
for these reasons, fpr evaluation is not pursued in this paper.
discussion
the duration hidden markov model proposed in this paper is a generic model for the oscillating structure of nucleosome and linker dnas in chromatin fiber.
the markov models can be replaced by any other models for the nucleosome and linker states.
the kernel method for linker length training is nonparametric and typically robust.
we showed in the simulation that updating the linker length distribution iteratively improves sensitivity and fdr in prediction if appropriate nucleosome and linker models are used.
in particular, the first iteration often achieves the most pronounced improvement.
in contrast inappropriate nucleosome and linker models could lead to the opposite outcome, as shown in the simulation studies .
in reality, the genomic dnas are complicated by their biological functions.
the models trained based on typical nucleosomes or linkers may not well fit some special genomic regions like repeated elements.
to avoid possible risks due to such complications, we trained the linker length distribution less greedily by using only one iteration in nupop.
limitations may still exist in the model training and assessment used in this study.
the mnase is known to have strong preference to cleave dinucleotides containing only a/t  <cit> .
consequently the mnase-mapped nucleosome sequences can be systematically biased in some regions.
this bias could undermine the prediction power because of the dampened signal in the trained nucleosome model.
the systematic bias may also exist in the well-defined nucleosomes, causing inaccuracy in sensitivity estimation.
a better map of nucleosomes is highly desirable for both purposes.
for species other than yeast, we currently lack high-quality genome-wide nucleosome sequence data  for model training and model validation.
the advantages of the re-scaling method shown using simulation in this paper need to be further assessed once such high-quality data becomes available.
moreover, the results from different methods in this paper were all based on the default settings.
the n-score method was originally trained based on a much smaller set of nucleosome and linker sequences.
a better training using a larger set could improve this method's predictions.
in addition, different settings in the n-score or mm/te methods can lead to different predictions, which we did not further investigate here.
finally, the software for mm/te method only provides the occupancy score.
different ways to call a predicted nucleosome based on the occupancy score might lead to different conclusions.
finally, we address the question of which subset of the available  <dig> reads data might best be used for training the nucleosome model.
in nupop, we trained the nucleosome model using selected non-redundant nucleosome reads of length within a short range , to retain strong high resolution nucleosome sequence signatures, e.g., the _ <dig> bp-periodic dinucleotide signals.
as comparisons, we trained two additional nucleosome models: one using the selected non-redundant reads of length 122- <dig> bp , and the other using all reads of length 122- <dig> bp.
the resulting models both contain the k-mer usage information that distinguishes nucleosomes from linkers , while the dinucleotide signals in these models are much weaker due to poor alignment of these reads.
furthermore, as the reads count at a nucleosome site is heavily biased by the g/c content due to mnase specificity and other effects in the experiment, the model trained from the redundant reads tends to be over-enriched in g/c.
when combined with the linker model from nupop, the two alternative nucleosome models yielded comparable sensitivity as nupop in predicting the approximate positioning of nucleosomes, assessed based on the  <dig>  well-defined nucleosomes used above.
this comparison, however, is not sensitive to spatial precision of the predictions.
therefore, we asked further, given that a nucleosome is predicted within ± <dig> of a true nucleosome, which model predicts the location more accurately?
to investigate this, we simulated genomic sequences using the nucleosome and linker models from nupop.
we compared the prediction from the three models and found that the true model with strong signals achieves much better prediction accuracy than the two alternative models.
for example,  <dig> % of the predictions from the true model were prefect , compared to  <dig> % and  <dig> % respectively from the other two models .
CONCLUSIONS
the dhmm model proposed in this paper is effective in characterizing the oscillating structure of nucleosome and linker dnas in chromatin fiber.
explicit modeling of linker length improves the prediction of nucleosome positioning regarding sensitivity.
the developed software tool nupop provides a user-friendly interface for predicting nucleosome occupancy and the most probable nucleosomes positioning map genome-wide.
availability and requirements
nupop software tools are freely available from http://nucleosome.stats.northwestern.edu.
the r package shall be made available through bioconductor http://www.bioconductor.org upon publication.
to run the nupop fortran stand-alone program, a fortran compiler is required.
for the nupop r package, an r version later than  <dig>  is required.
authors' contributions
lxi did all the data analyses.
jpw, jw, and lxi wrote the paper.
jpw and jw directed the research.
yfm conducted all lab work for data generation and validation.
jpw and lxi developed the nupop software tools.
lxia and jf implemented the nupop web-server.
all authors read and approved the final manuscript.

 most eukaryotic genomic dna is wrapped in nucleosomes, which occlude and strongly distort the wrapped dna. accumulating evidence shows that the dna sequence itself is highly predictive of nucleosome positioning in vivo, and that nucleosome positioning is closely related to chromosome functions. a fast software tool for predicting nucleosome positioning is highly desirable. several statistical methods for nucleosome positioning prediction have been proposed in the literature. in a method was proposed based on cross-correlation with a nucleosome dna sequence signature. in a markov model was used together with consideration of steric exclusion and thermodynamic equilibrium. in, a support vector machine was trained based on the differential k-mer usage between nucleosome and linker dnas. in, the authors proposed an n-score model to discriminate nucleosome and linker dnas using wavelet energies as covariates in a logistic regression model. in, a web-interface called "nuscore" for estimation of the affinity of histone core to dna and prediction of nucleosome positioning was developed based on the dna deformation energy score. in, the model from was improved by incorporation of differential k-mer usage tracts, which are strongly disfavored by nucleosomes). this model can be further improved by accounting for nucleosome-nucleosome interaction. while the nucleosomal features are universal, eukaryotic genomes vary in nucleosomal repeat length and base composition. the nucleosomal repeat length is dictated by the length distribution of linker dnas that separate neighboring nucleosomes, and it determines the overall nucleosome density in the chromatin fiber. the contribution of this paper is a duration hidden markov model and a software tool called nupop for genome-wide nucleosome positioning prediction. we show that incorporation of linker length information can achieve better sensitivity in prediction. in addition, we propose a re-scaling method to adjust for base composition variation when using yeast models to make predictions for other species. a relatively superior performance of this approach is established by comparing it with other existing tools. methods model the hidden markov model has been known for decades. an excellent and famous tutorial is rabiner's paper, in which the model, algorithms, and applications were carefully and thoroughly reviewed. a conventional hmm implicitly assumes a geometric duration distribution for each state, which can be wrong in real applications. modeling the explicit duration of each state can improve the prediction of hmms. we model each chromosomal dna sequence with a duration hidden markov model of two oscillating states: nucleosome and linker, where the nucleosome state has a fixed length of bp, and the linker state has a variable length. we assume that at the end of each state, the chain must transit to the other state; additionally, a complete chromatin sequence must start with and end in a linker state. we trained a 4th order time-dependent markov chain for the the n state, and a homogeneous 4th order markov chain for the l state to distinguish the k-mer usage preferences for k up to between the nucleosome and linker states as shown in other methods, e.g.,. let e = e..., e be a nucleosome dna sequence. let pn be the probability of observing e as a nucleosome, computed as the product of probabilities for both watson and crick strands under the 4th order markov chain model. we assume that the linker dna length of a given species has an unknown distribution fl defined for k =..., τl. an observed linker dna sequence e = e..., ek carries two pieces of information, the length is k bp, and given which, the emitted letters are e..., ek. let gl denote the homogeneous markov chain model for the linker dna. then observing e as a linker dna has probability suppose x = x..., xn is a genomic dna sequence of length n, where xi = a/c/g/t. let z = z..., zn be the corresponding hidden state path, where zi = if xi is covered by a nucleosome state, and otherwise. suppose that the path z = z..., zn partitions x into k consecutive nucleosome or linker state blocks, in which the nucleosome blocks have a uniform length of bp, whereas the length of linker blocks may vary. we denote these blocks as y = y+,..., yb, and their state identification as s = s..., sb, where si = if yi is nucleosome state, and otherwise. the probability of observing is given by where π and πe stand for the probabilities that the chain initializes and ends with state s and sk respectively, and i is an indicator function. since we assume that a complete chromatin sequence must start with and end in a linker state, π = πe = we define the nucleosome occupancy at a specific position i, denoted oi, as the posterior probability that zi = i.e., we also define the histone binding affinity score at position i as the log likelihood ratio for the region xi-... xi,..., xi+ to be a nucleosome vs. a linker, i.e., given the models pn, gl and fl, the optimal path z can be found by the standard viterbi algorithm, and the nucleosome occupancy score can be estimated using forward and backward algorithms. data and model training we utilized the yeast nucleosome dna reads from pyrosequencing published in for model training and assessment. among reads that each were mapped to a unique region of the yeast genome, we first selected reads of length between and bp. if multiple such reads existed for the same nucleosome, we selected the one with the highest blast score. the resulting non-redundant set of nucleosome sequences were center aligned to train the nucleosome model pn. the 4th order time dependent markov chain can be defined by the base composition at the first position qn, and the transitional probabilities qn, qn, qn, qn, for k =..., xi = a/c/g/t, i =..., where the subscript k, i index the positions within a nucleosome. these probabilities are trained using the corresponding observed fractions or conditional fractions based on the center alignment, with a three bp moving average. we further identified reads-free regions of length 7- bp to train the linker state model gl. the 4th order homogeneous markov model for the linker dnas can be completely defined by the stationary base composition ql, and the transition probabilities ql, ql, ql, ql. by "homogeneous", we mean that these probabilities are all constants as functions of i. these probabilities were trained using their observed values as in the nucleosome model. for example, ql was trained by calculating the fraction of occurrences of transitions from any four letters to the fifth letter in the selected putative linker dna sequences. our initial nucleosome/linker model was trained using the yeast data. a complication arises when predicting nucleosomes for other species because organisms may differ significantly in their dna base composition. we propose to scale up or down the probabilities in the markov models by a factor determined by the difference of the base composition between the current species and yeast. for example, in c. elegans, the fraction of a plus t bases is compared to in yeast. for a specific transition probability qn at any specific nucleosomal position defined for yeast, we scaled it up as qn × / for the corresponding transition probability at that given position for c. elegans. likewise the transition probabilities for g/c will be scaled down by a factor of /. all the re-scaled probabilities are then normalized. the same re-scaling applies to the linker model. we shall use simulations below to show that re-scaling improves prediction regarding sensitivity and false discovery rate. using the trained nucleosome model and linker model, we further train the linker dna length distribution as follows. we assume that the linker dnas in any given species or cell type have a maximum length τl = bp. this algorithm contains the following steps: initialize the algorithm with a uniform distribution for fl for k =... τl where τl is the maximum allowable linker length. use the forward and backward algorithm to obtain the posterior expectation of fl for each k. for a sequence x = x..., xn, let nk be the number of linker dnas of length k. then for k =..., τl. for k =..., τl. update the empirical linker length distribution from step using a kernel smoothing method as follows: where k is the standard gaussian kernel and h is the bandwidth parameter optimally chosen by the leave-one cross-validation method as in. use the updated linker length distribution from step to compute the nucleosome occupancy and optimal positioning. compared to viterbi training, using the posterior expectation obtained in eq. combined with the kernel method in eq. performs overwhelmingly better in minimizing the summed square errors. in the developed software tool nupop, we have trained the linker dna length distributions for different species including human, mouse, rat, zebrafish, d. melanogaster, c. elegans, s. cerevisiae, c. albicans, s. pombe; a. thaliana and maize. the linker dna length distribution for each species has been trained by scanning the corresponding genome sequences based on τ l = we found that the re-scaled nucleosome and linker profiles, together with the trained linker length distribution, not only roughly recover the genome-wide base compositions, but also the dinucleotide frequencies for different species. the frequency of each single or di-nucleotide in simulated genomes typically differs by ≤ 1% from that observed in the corresponding real genomes. as different cell types from the same organism can exhibit quite different linker dna length distributions, a useful future refinement would utilize high quality nucleosome maps for the given cell type, when such data become available. software tools we have developed a software tool called nupop, implemented in three different formats: an r package tested for windows xp, linux and mac os x; a stand-alone fortran program; and an nupop web server, all available from http://nucleosome.stats.northwestern.edu. the r package is built upon the fortran program. it provides additional handy functions to visualize the resulting viterbi and nucleosome occupancy predictions. both the r package and fortran program can handle a genomic sequence of any length with a ram demand < m bytes. the predicted results are stored locally in the working directory. the web server provides an interface through which the user can submit their own sequence up to k bp in length for fast online prediction. when making a prediction, the user is required to specify which species the genomic sequence is from. if the species is not on the list, nupop will calculate the base composition of the input dna sequence and then choose the nucleosome and linker models from a species that has the most similar base composition. an alternative model with a 1st order time-dependent markov chain for the nucleosome state and a homogeneous 1st order markov model for the linker state, trained in the same way, is also implemented in nupop as an option.  updating fl improves prediction updating the linker length distribution not only helps recover the true nucleosome density, but also improves prediction. we demonstrated this by simulation as follows. we simulated genomic sequences with the 4th and 1st order yeast models respectively, each containing nucleosomes and linkers. the linker dna length was simulated from a normal distribution and a gamma distribution. if a nucleosome is predicted within ± bp of a true nucleosome, we define it as a correct positive prediction. the rate of the correct positive prediction, referred to as sensitivity, is defined as the percentage of the nucleosomes that are correctly predicted. in addition, we include the false discovery rate, defined as the fraction of the predicted nucleosomes that reside > ± bp away from any true nucleosomes, as the second measure for model performance. in analogy to statistical hypothesis testing, the sensitivity measures the power of prediction, while the fdr measures the fraction of type i errors in the positive claims. the results are presented in table and in both cases, the linker length was initialized as a uniform distribution with τl = compared to the gamma distribution, the normal distribution is relatively flatter. thereby updating the linker length distribution did not significantly change the total number of predicted nucleosomes. the sensitivity increased on average by ~ 4-5% and the fdr dropped by ~ 5% after one update. further updating continued to improve the prediction until it stabilized after four iterations. in contrast, the gamma model is much more skewed. the uniform linker length distribution resulted in an under-estimation of the total number of nucleosomes. by four updatings, the sensitivity increased by 8%, while fdr remained relatively more stable. total predictions, sensitivity, and false discovery rate are the averages based on repeated simulations. for each simulation a genomic sequence consisting of nucleosomes and linkers were simulated using the 1st and 4th order yeast models. the linker length distribution was initialized as uniform on..., and was iteratively updated in the dhmm. results are shown after and updates. total predictions, sensitivity, and false discovery rate are the averages based on repeated simulations. for each simulation a genomic sequence consisting of nucleosomes and linkers were simulated using the 1st and 4th order yeast models. the linker length distribution was initialized as uniform on..., and was iteratively updated in the dhmm. results are shown after and updates. we also observe that under the same setting and condition, the fourth order model performs slightly but uniformly better than the first order model in both sensitivity and fdr. this is given that the true models are known and we scan the sequence using the true models. in theory, the first order markov chain model is nested in the fourth order model. therefore if the true model is the first order, a well trained fourth order model will have the same prediction power as the first order model, but not vice versa. since training a higher order markov chain model requires more data, inadequate training can undermine the prediction power. re-scaling vs. not re-scaling to illustrate the advantages of re-scaling, we re-scaled the yeast profiles according to the base composition of the maize genome. using the re-scaled profiles we simulated genomic sequences that each contain nucleosomes and linkers. the linker dna length followed the same two distributions as in table and we compare the prediction results from the scaled and non-scaled yeast profiles in table and we found using the re-scaled profile yields a lower fdr than using the yeast profile. in addition, updating the linker length under the correct profile consistently improves the sensitivity and fdr until prediction stabilizes. in contrast, while using the yeast profile to scan the simulated maize-like genome, the prediction drastically deteriorates as the linker length updating proceeds. the same simulation was repeated on other species including human and c. elegans, where the base composition is similar to yeast. unsurprisingly, the results from the scaled yeast profile were still better than those from the original yeast profile in terms of both sensitivity and fdr, while the difference is much smaller than for the maize case. total predictions, sensitivity, and false discovery rate are the averages based on repeated simulations. for each simulation a maize-like genomic sequence consisting of nucleosomes and linkers were simulated using the re-scaled 1st and 4th order yeast models. each sequence was scanned using the true models and the yeast models with an initial uniform linker length distribution on..., the results after updates of linker length distribution are compared. total predictions, sensitivity, and false discovery rate are the averages based on repeated simulations. for each simulation a maize-like genomic sequence consisting of nucleosomes and linkers were simulated using the re-scaled 1st and 4th order yeast models. each sequence was scanned using the true models and the yeast models with an initial uniform linker length distribution on..., the results after updates of linker length distribution are compared. nupop vs. other software tools we briefly assess the prediction performance of nupop by comparing it with two existing methods: the n-score method of and the markov model/thermodynamic equilibrium method of. as the exact genome-wide nucleosome positioning map is unknown, we utilize the high-throughput sequence reads to identify well-defined nucleosomes. we first selected sequences of length between 130- bp and constructed a reads-based occupancy map. the reads-occupancy score at a specific position is defined as the number of reads that covered this position. then we calculated the moving average of this reads occupancy score using a bp window. a sharp peak in the average occupancy curve indicates a nucleosome with well-defined positioning. considering that the average linker dna length is bp in yeast, we quantified the sharpness of the peak by calculating the slope from the peak point to the up-/down-stream bp point on the average occupancy curve. we set one condition for the peak to be selected as the center of a nucleosome to be that the absolute value of slope from either side should be >. secondly, we required that the peak height itself must be ≥, i.e., a well-defined nucleosome must be testified by at least two well-overlapped reads. we chose the threshold value as instead of because the overlap of the two reads can be less than bp, resulting in a peak on the moving average curve slightly lower than. with these criteria, a total of well-defined nucleosomes are selected from the chromosomes of yeast. a snapshot of a region with many selected well-defined nucleosomes is presented in figure figure provides a snapshot of nucleosome occupancy predicted by nupop together with the reads-occupancy. we first assess the sensitivity of predictions from nupop using the well-defined control set. if there is a predicted nucleosome within ±k bp of any well-defined nucleosomes, we count this as one correct prediction. we varied k from... to investigate the sensitivity behavior at different precision thresholds. the n-score model predicted nucleosomes. the current software tool for mm/te method does not provide viterbi predictions, but only the nucleosome occupancy scores. therefore we calculated the moving average of the occupancy score using a bp moving window. the resulting peaks were treated as the centers of predicted nucleosomes. if two peaks reside within bp, we discarded the one with smaller moving average of occupancy score. this procedure identified predicted nucleosomes. likewise, using the occupancy scores from nupop, we identified and predicted nucleosomes under the 4th and 1st order models respectively. in figure 3a, we compare the sensitivity estimates from the 4th order model of nupop with the n-score and mm/te methods at different threshold values of prediction accuracy. as the sensitivity tends to increase with an increase in the total predictions, we further selected and best predictions from nupop model to compare with the n-score method and in figure 3b and 3c respectively. the sensitivity results suggest that the predictions from nupop outperforms the other two methods in two senses. firstly, the sensitivity from nupop is - % higher than the other two methods at different threshold values. secondly, while controlling the total predictions to be the same, nupop has ~ - % higher sensitivity than the mm/te method, when the precision threshold is ≤ ± as the precision threshold gets less stringent, the difference attenuates and eventually vanishes. the contrast between nupop and the n-score method is even larger as shown in figure 3c. as a further comparison, we computed the predictions of the dhmm method under a uniform linker length distribution defined on..., this method predicted nucleosomes under the 4th order models, achieving a sensitivity - % higher than the mm/te method. when controlling the total predictions to be the same as mm/te method or n-score method, the resulting sensitivity curve almost perfectly overlaps with that from nupop. therefore we omitted these results from figure 3b and 3c. one could further attempt to evaluate the false positive rate, measuring the fraction of linker regions that were falsely classified as nucleosome regions. this task requires well-defined linker regions. a problem, however, is that the average length of linker dnas in yeast is smaller than the dispersion in lengths of the nucleosome dnas as isolated biochemically. thus existing nucleosome maps lack the precision needed to define such short linker dnas. moreover, various sampling biases such as the dna sequence preferences of the micrococcal nuclease used to liberate nucleosomes biochemically could yield longer genomic regions that are free of recovered nucleosome dna reads even if they are actually nucleosome occupied. attempts to evaluate the fpr given these problems in the data could result in misleading conclusions. for these reasons, fpr evaluation is not pursued in this paper. discussion the duration hidden markov model proposed in this paper is a generic model for the oscillating structure of nucleosome and linker dnas in chromatin fiber. the markov models can be replaced by any other models for the nucleosome and linker states. the kernel method for linker length training is nonparametric and typically robust. we showed in the simulation that updating the linker length distribution iteratively improves sensitivity and fdr in prediction if appropriate nucleosome and linker models are used. in particular, the first iteration often achieves the most pronounced improvement. in contrast inappropriate nucleosome and linker models could lead to the opposite outcome, as shown in the simulation studies. in reality, the genomic dnas are complicated by their biological functions. the models trained based on typical nucleosomes or linkers may not well fit some special genomic regions like repeated elements. to avoid possible risks due to such complications, we trained the linker length distribution less greedily by using only one iteration in nupop. limitations may still exist in the model training and assessment used in this study. the mnase is known to have strong preference to cleave dinucleotides containing only a/t. consequently the mnase-mapped nucleosome sequences can be systematically biased in some regions. this bias could undermine the prediction power because of the dampened signal in the trained nucleosome model. the systematic bias may also exist in the well-defined nucleosomes, causing inaccuracy in sensitivity estimation. a better map of nucleosomes is highly desirable for both purposes. for species other than yeast, we currently lack high-quality genome-wide nucleosome sequence data for model training and model validation. the advantages of the re-scaling method shown using simulation in this paper need to be further assessed once such high-quality data becomes available. moreover, the results from different methods in this paper were all based on the default settings. the n-score method was originally trained based on a much smaller set of nucleosome and linker sequences. a better training using a larger set could improve this method's predictions. in addition, different settings in the n-score or mm/te methods can lead to different predictions, which we did not further investigate here. finally, the software for mm/te method only provides the occupancy score. different ways to call a predicted nucleosome based on the occupancy score might lead to different conclusions. finally, we address the question of which subset of the available reads data might best be used for training the nucleosome model. in nupop, we trained the nucleosome model using selected non-redundant nucleosome reads of length within a short range, to retain strong high resolution nucleosome sequence signatures, e.g., the _ bp-periodic dinucleotide signals. as comparisons, we trained two additional nucleosome models: one using the selected non-redundant reads of length 122- bp, and the other using all reads of length 122- bp. the resulting models both contain the k-mer usage information that distinguishes nucleosomes from linkers, while the dinucleotide signals in these models are much weaker due to poor alignment of these reads. furthermore, as the reads count at a nucleosome site is heavily biased by the g/c content due to mnase specificity and other effects in the experiment, the model trained from the redundant reads tends to be over-enriched in g/c. when combined with the linker model from nupop, the two alternative nucleosome models yielded comparable sensitivity as nupop in predicting the approximate positioning of nucleosomes, assessed based on the well-defined nucleosomes used above. this comparison, however, is not sensitive to spatial precision of the predictions. therefore, we asked further, given that a nucleosome is predicted within ± of a true nucleosome, which model predicts the location more accurately? to investigate this, we simulated genomic sequences using the nucleosome and linker models from nupop. we compared the prediction from the three models and found that the true model with strong signals achieves much better prediction accuracy than the two alternative models. for example, % of the predictions from the true model were prefect, compared to % and % respectively from the other two models.  the dhmm model proposed in this paper is effective in characterizing the oscillating structure of nucleosome and linker dnas in chromatin fiber. explicit modeling of linker length improves the prediction of nucleosome positioning regarding sensitivity. the developed software tool nupop provides a user-friendly interface for predicting nucleosome occupancy and the most probable nucleosomes positioning map genome-wide. availability and requirements nupop software tools are freely available from http://nucleosome.stats.northwestern.edu. the r package shall be made available through bioconductor http://www.bioconductor.org upon publication. to run the nupop fortran stand-alone program, a fortran compiler is required. for the nupop r package, an r version later than is required. authors' contributions lxi did all the data analyses. jpw, jw, and lxi wrote the paper. jpw and jw directed the research. yfm conducted all lab work for data generation and validation. jpw and lxi developed the nupop software tools. lxia and jf implemented the nupop web-server. all authors read and approved the final manuscript. 
nupop provides a user-friendly software tool for predicting the nucleosome occupancy and the most probable nucleosome positioning map for genomic sequences of any size. simulation studies show that modeling the linker length distribution and utilizing a base composition re-scaling method both improve the prediction of nucleosome positioning regarding sensitivity and false discovery rate. therefore a fast software tool for predicting nucleosome positioning can help understanding how a genome's nucleosome organization may facilitate genome function. the nucleosome and linker models trained from yeast data are re-scaled when making predictions for other species to adjust for differences in base composition. we present a duration hidden markov model for nucleosome positioning prediction by explicitly modeling the linker dna length. increasing evidence has shown that genomic dna sequence itself is highly predictive of nucleosome positioning genome-wide. a software tool named nupop is developed in three formats for free download.
 the nucleosome is the fundamental packing unit of dnas in eukaryotic cells. its detailed positioning on the genome is closely related to chromosome functions. increasing evidence has shown that genomic dna sequence itself is highly predictive of nucleosome positioning genome-wide. therefore a fast software tool for predicting nucleosome positioning can help understanding how a genome's nucleosome organization may facilitate genome function.  we present a duration hidden markov model for nucleosome positioning prediction by explicitly modeling the linker dna length. the nucleosome and linker models trained from yeast data are re-scaled when making predictions for other species to adjust for differences in base composition. a software tool named nupop is developed in three formats for free download.  simulation studies show that modeling the linker length distribution and utilizing a base composition re-scaling method both improve the prediction of nucleosome positioning regarding sensitivity and false discovery rate. nupop provides a user-friendly software tool for predicting the nucleosome occupancy and the most probable nucleosome positioning map for genomic sequences of any size. when compared with two existing methods, nupop shows improved performance in sensitivity. 