Predicting nucleosome positioning using a duration Hidden Markov Model

Model

The hidden Markov model (HMM) has been known for decades. An excellent and famous tutorial is Rabiner's 1989 paper [14], 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 (e.g., [14–17]). We model each chromosomal DNA sequence with a duration hidden Markov model (dHMM) of two oscillating states: nucleosome (N) and linker (L), where the nucleosome state has a fixed length of 147 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 5 between the nucleosome and linker states as shown in other methods, e.g., [3, 6] (see below for details).

Let e = e₁, ..., e₁₄₇ be a nucleosome DNA sequence. Let P _N 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 F _L (k) defined for k = 1, ..., τ _L (the maximum linker length we allow). An observed linker DNA sequence e = e₁, ..., e _k carries two pieces of information, the length is k bp, and given which, the emitted letters are e₁, ..., e _k . Let G _L (e|k) denote the homogeneous Markov chain model for the linker DNA (again including both strands). Then observing e as a linker DNA has probability

Suppose x = x₁, ..., x _n is a genomic DNA sequence of length n, where x _i = A/C/G/T. Let z = z₁, ..., z _n be the corresponding hidden state path, where z _i = 1 if x _i is covered by a nucleosome state, and 0 otherwise. Suppose that the path z = z₁, ..., z _n partitions x into k consecutive nucleosome or linker state blocks, in which the nucleosome blocks have a uniform length of 147 bp, whereas the length of linker blocks may vary. We denote these blocks as y = y+, ..., y _B , and their state identification as s = s₁, ..., s _B , where s _i = 1 if y _i is nucleosome state, and 0 otherwise. The probability of observing (x, z) is given by

where π₀(s₁) and πe(s _B ) stand for the probabilities that the chain initializes and ends with state s₁ and s _k respectively, and I is an indicator function. Since we assume that a complete chromatin sequence must start with and end in a linker state, π₀(s₁ = 0) = π _e (s _B = 0) = 1. We define the nucleosome occupancy at a specific position i, denoted o _i , as the posterior probability that z _i = 1, i.e.,

We also define the histone binding affinity score at position i as the log likelihood ratio for the region x _i -73, ... x _i , ..., x _i +73 to be a nucleosome vs. a linker, i.e.,

Given the models P _N , G _L and F _L , 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 503,264 yeast nucleosome DNA reads from 454 pyrosequencing published in [6] for model training and assessment. Among 371,914 reads that each were mapped to a unique region of the yeast genome, we first selected reads of length between 146 and 149 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 18,547 nucleosome sequences were center aligned to train the nucleosome model P _N . The 4th order time dependent Markov chain can be defined by the base composition at the first position qN(x₁), and the transitional probabilities qN(x₂|x₁), qN(x₃|x₁, x₂), qN(x₄|x₁, x₂, x₃), qN(x _k |x_k-4, x _k -3, x_k-2, x_k-1), for k = 5, ..., 147, x _i = A/C/G/T, i = 1, ..., 147, 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 (as explained in [1, 18]). We further identified 8,090 reads-free regions of length 7-500 bp to train the linker state model G _L . The 4th order homogeneous Markov model for the linker DNAs can be completely defined by the stationary base composition qL(xi), and the transition probabilities qL(x _i |x_i-1), qL(x _i |x_i-1, x_i-2), qL(x _i |x_i-1, x_i-2, x_i-3), qL(x _i |x_i-1, x_i-2, x_i-3, x_i-4). 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(x _i |x_i-1, x_i-2, x_i-3, x_i-4) 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 0.645 compared to 0.617 in yeast. For a specific transition probability qN(A|....) at any specific nucleosomal position defined for yeast, we scaled it up as qN(A|....) × 0.645/0.617 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 0.355/0.383.

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 (P _N ) and linker model (G _L ), 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 = 500 bp.

for k = 1, ..., τ _L .

(1)

Compared to Viterbi training (i.e., using linker length predicted from the Viterbi algorithm), using the posterior expectation obtained in Eq. (1) combined with the kernel method in Eq. (2) performs overwhelmingly better in minimizing the summed square errors (unpublished work [17]). In the developed software tool NuPoP, we have trained the linker DNA length distributions for 11 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 (F _L ) for each species has been trained by scanning the corresponding genome sequences based on τ _L = 500. 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 (results not shown). As different cell types from the same organism (with the same genome) can exhibit quite different linker DNA length distributions [13], 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 (most probable nucleosome position map) and nucleosome occupancy predictions. Both the R package and Fortran program can handle a genomic sequence of any length with a RAM demand <400 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 500 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 F _L improves prediction

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

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 [4] (results kindly provided by Dr. G. Yuan, personal communication) and the Markov model/thermodynamic equilibrium method of [7] (to be called MM/TE method below). As the exact genome-wide nucleosome positioning map is unknown, we utilize the 371,914 454 high-throughput sequence reads to identify well-defined nucleosomes. We first selected sequences of length between 130-160 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 147 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 20 bp in yeast [13], we quantified the sharpness of the peak by calculating the slope from the peak point to the up-/down-stream 20 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 > 0.01. Secondly, we required that the peak height itself must be ≥ 1.9, i.e., a well-defined nucleosome must be testified by at least two well-overlapped reads. We chose the threshold value as 1.9 instead of 2.0 because the overlap of the two reads can be less than 147 bp, resulting in a peak on the moving average curve slightly lower than 2.0. With these criteria, a total of 20,471 well-defined nucleosomes are selected from the 16 chromosomes of yeast. A snapshot of a region with many selected well-defined nucleosomes is presented in Figure 1. Figure 2 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 (center to center), we count this as one correct prediction. We varied k from 5, 10, ... 70, 73 to investigate the sensitivity behavior at different precision thresholds. The N-score model predicted 48,394 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 147 bp moving window. The resulting peaks were treated as the centers of predicted nucleosomes. If two peaks reside within 127 bp, we discarded the one with smaller moving average of occupancy score. This procedure identified 43,979 predicted nucleosomes (if we had required two nucleosomes to be 147 bp away, even fewer nucleosomes would have been identified). Likewise, using the occupancy scores from NuPoP, we identified 52,327 and 51,380 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 (the fourth order model from NuPoP performs better, but very slightly, than the first order model. Hence the latter was omitted in Figure 3 for better presentation). As the sensitivity tends to increase with an increase in the total predictions, we further selected 43,979 and 48,394 best predictions from NuPoP model (here "best" is in the sense of the largest sum of occupancies over 147 bp) 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 4.9-8.3% higher than the other two methods at different threshold values (Figure 3a). Secondly, while controlling the total predictions to be the same, NuPoP has ~ 3.2-5.3% higher sensitivity than the MM/TE method (Figure 3b), when the precision threshold is ≤ ±35. 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 1, 2,..., 500. This method predicted 48,334 nucleosomes under the 4th order models, achieving a sensitivity 3.0-6.6% higher than the MM/TE method (Figure 3a). 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 (FPR), measuring the fraction of linker regions that were falsely classified as nucleosome regions (or similarly the false discovery rate, FDR). This task requires well-defined linker regions. A problem, however, is that the average length of linker DNAs in yeast (20 bp; [13]) is smaller than the dispersion in lengths of the nucleosome DNAs as isolated biochemically (which is often 30-50 bp full width at half maximum, notwithstanding that the nucleosome as defined crystallographically has precisely 147 bp of DNA). 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 (which preferentially cleaves A/T rich regions) could yield longer genomic regions that are free of recovered nucleosome DNA reads even if they are actually nucleosome occupied [20]. 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.