New decoding algorithms for Hidden Markov Models using distance measures on labellings

HMM definitions

A hidden Markov Model (HMM) is a tuple M = (A, E, Σ, x₀): A is the m × m transition matrix where a _ij gives the probability of transition from state i to state j; E is an m × |Σ| emission probability matrix where e _kσ is the probability of emitting symbol σ ∈ Σ in state k, and x₀ is the start state of the model. A path in an HMM is a sequence of states x₀, x₁, ..., x _n ; in step i of the execution of the model, we transition from x_i-1to a new state x _i , and emit symbol y _i according to the distribution of row x _i of the matrix E. The HMM defines a probability measure over paths and sequences: the joint probability of sequence y = y₁, ..., y _n and path x = x₀, x₁, ..., x _n is .

In a labelled HMM, we add a labelling function ℓ, which assigns to each state of the model (from 1 to m) a label, which typically corresponds to a sequence feature. For a path x, let λ = λ₀, λ₁, ..., λ _n = ℓ(x₀), ℓ(x₁), ...,ℓ(x _n ) be its labelling. Many states may share a label, so many paths may also share a labelling.

In HMM decoding, we are interested in labellings, which assign a feature to each position of a sequence. Often, a labelling for a sequence will have many consecutive positions with the same label. Given a labelling λ = λ₀, f₁, f₁, ..., f₁, f₂, f₂, ..., f₂, ..., f _k , ..., f _k , which consists of λ₀ followed by a number of positions labelled f₁, then a number of positions labelled f₂ (which is different from f₁), and so on, we define its footprint to be the sequence f = f₁, ..., f _k ; this corresponds to the overall labelling of the sequence, but with the feature boundaries left entirely flexible. For a sequence of length n, its footprint may be much shorter than n in length, assuming it has long features.

For each label z, let L _z be the set of states with label z; let L be the size of the largest L _z .

Distance measures for labellings

Here, we consider two different types of distance measures for labellings of the same sequence. The appropriate distance measure depends on the domain: our first type is more suited to short footprints corresponding to long features, while the second, based on Hamming distance, is more suited to long footprints.

Consider two different labellings λ and λ' for the same sequence y. We will define their border shift distance d(λ, λ') as a function of the boundaries of their features.

First, if λ and λ' have different footprints, they are incompatible; we will set their distance to ∞.

Otherwise, let f = f₁, ..., f _k be their common footprint, and let b _i (λ) be the position in λ corresponding to the first emission from the interval labelled f_i+1. The border shift distance between λ and λ' is then max_{i = 1... k-1}|b _i (λ) - b _i (λ')|: the maximum boundary shift between the corresponding features in the two labellings. Figure 1 shows how we compute this distance.

We also define the border shift sum distance, which is infinite for unequal footprints, and otherwise.

A quite different measure is based on counting mismatches between the labellings. In the generalized Hamming distance, we first assume we have a non-negative integer-valued distance matrix H that gives the distance between any two labels i and j. (In the easiest case, H _ij = 0 if i = j and 1 otherwise; this gives the Hamming distance.) Then, is the generalized Hamming distance between the two labellings.

For any metric d on labellings, we define the ball of radius r around a centre λ: B _d (λ, r) = {λ' : d(λ, λ') ≤ r}.

With this in mind, we state our goal. Given an HMM M, a sequence y, a distance metric d, and a distance tolerance r, we seek the labelling λ* such that Pr[B _d (λ*, r)|y, M] is maximized. This labelling has, according to M's probability measure, the highest probability of having the true labelling for the sequence be at most distance r away, according to our distance measure.

Traditional algorithms

Algorithms typically used for decoding HMMs optimize quite different goals from ours. The most common are the Viterbi algorithm [3], which maximizes the probability of the state path used to generate a sequence, or the posterior algorithm, which maximizes the expectation of the number of positions in the sequence assigned to the correct state.

Both of these algorithms experience serious problems. The Viterbi path may itself have extremely tiny probability, and for long sequences, is highly unlikely to be correct overall; also, since this algorithm maximizes a global criterion, its prediction may be locally very sloppy. The posterior algorithm gives credible guesses for each position by itself, but does not identify intervals of the sequence as having common features; indeed, adjacent positions may be labelled with states that cannot co-exist in paths in the model with nonzero probability.

One might attempt to use generalized HMMs in this context, where the length distribution of a state is no longer geometric. However, this does nothing to allow the boundary between labels to be flexible, and may not help when multiple states share the same label.

One approach to hybridizing these methods is to compute the path of nonzero probability that maximizes the posterior decoding objective [4], or that maximizes the geometric mean of the probability of assigning correct states to positions [5]. These algorithms may improve results, but depend on the topology of the model being correct, and also may still give paths that are globally quite poor. There are variants of them that also work with labellings, not paths.

A tempting alternative might be to compute the most probable labelling of a sequence instead of the Viterbi path; this corresponds to our objective, but for balls of radius zero. This objective is NP-hard to optimize [6] in general, though Brejová et al., in addition to proving the NP-hardness of the overall problem, do give an optimization algorithm for some special cases. Krogh [1] proposed the 1-best algorithm, a heuristic to find a highly probable labelling, which is used in many applications, including transmembrane protein topology prediction [2]. However, even the most probable labelling may have very small conditional probability.

Several authors have also considered sets of paths rather than individual paths. For example, we could sample k paths from an HMM to predict alternatively-spliced exons [7, 8]. In recent work, Brown and Golod [9] proposed finding the k most probable paths through an HMM, and applied this idea to transmembrane protein topology prediction. Our methods consider an exponential number of implicitly represented labellings, in contrast to these methods.

Computing the probability of a labelling and of a footprint

We begin with the simple problems of computing the probability of a labelling λ or of a footprint f. For both of these problems, our algorithm is a variant of the traditional forward algorithm for HMMs [3], which computes the probability that an HMM generates a sequence y in O(nm²) time.

To compute the probability of all paths compatible with the labelling, we follow the usual forward algorithm, but at position i in the sequence, we restrict to allowing only transitions in the model between states labelled λ_i-1and λ _i ; other transitions are given probability zero. This algorithm takes less time than the forward algorithm to run; if the largest number of edges between any two (possibly identical) labels λ _i and λ _j is q, then this algorithm takes O(qn) time, which is O(L²n).

To compute the probability of all paths compatible with the footprint f = f₁, ..., f _k , we follow the usual forward algorithm, but on a variant HMM. We will create k groups of states G₁, ..., G _k , each corresponding to a label set , for each position f _i in the footprint: each state in is represented in G _i . States in the new model have exactly the same emission probabilities as in M, but we may only make transitions from states in G _i to those in G _i or in G_i+1, with the same transition probabilities as in M . (To make a proper HMM, we can create a "dump" state for paths in M that do not respect the footprint we seek, or not bother.) See Figure 2.

We require that the initial transition is from x₀ to a state in G₁, and we compute only the sum of the forward probabilities ending in states in G _k that have made all of the required label transitions in the footprint f. The new model has O(Lk) states and each state has degree at most 2L, so computing the probability of a footprint can be done using the forward algorithm in O(L²kn) time. Note that this could be quadratic in n, as k, the length of the footprint, is bounded above by n.

Computing the probability of a ball of labellings

Computing the probability of B(λ, r), the ball of labellings around λ with radius r can be done efficiently for the border shift distance. As before, we produce a group of states G _i for each position f _i in the footprint f = f₁ ... f _k of labelling λ . Each state in G _i has nonzero transition probability only to states in G _i or G_i+1. However, if b _i (λ) is the position of the i-th border in the labelling λ, then we allow states in G _i to appear only in positions from b_i-1- r to b _i + r - 1. For each position j, we can identify the groups that can be active at that position in this way. Computing the ball probability is easily accomplished using the pattern for computing the probability of a labelling: we simply set the transition probabilities for transitions to a new label that are found outside of the correct interval to zero and forbid transitions within a label after its valid interval ends. We then run the forward algorithm on the updated HMM, which we call M'. M' has exactly the same states as for computing a footprint, meaning that the runtime to compute the probability of the ball B(λ, r) is O(L²kn). However, we can reduce this runtime. If t _j is the number of groups of states G _i that could be active at position j, then allowing the boundaries to shift by r makes t _j go up by 1 at 2r positions. The asymptotic runtime is q∑t _j , which is O(qn + qrk). Again, since q ≤ L², this is O(nL² + rkL²), which may be well less than the runtime of the forward algorithm on the original model M, depending on r and k.

For the border shift sum distance, we must remember the total amount by which we have already shifted the first j - 1 borders before we handle the next border. As such, we can imagine the same sort of HMM as for the maximum border shift, but with an additional parameter in each state to keep track of this total mismatch so far accumulated. The model winds up with r times as many states this way, but the set of active groups at each position, and the degree of each node in the new HMM, stay the same. The overall runtime is r times slower than for the maximum border shift distance, or O(qnr + qr²k). We note also that the values of r chosen in practical applications are likely to be much larger for this distance measure than for the maximum border shift measure; as such, we expect that these algorithms are not likely to be as useful in practice.

Finally, for the generalized Hamming distance, with label distance matrix H, we can again hybridize the forward algorithm. Consider again our reference labelling λ, and create a new HMM with (r + 1)m states, where each state's identifier is (j, d). The semantic meaning of being in state (j, d) at position i of the sequence is that the HMM M is in state j and we have accumulated d distance from λ in the first i positions of the sequence; we create such states for all d in the range from 0 to r. If there is a transition from j to k in the HMM M, and if H_{L(λ{i=1}), L(k)}= h, then we draw an edge in our new HMM from (j, d) to (k, d + h) for all d such that d + h ≤ r; if d + h is more than r, then we do not allow the transition. The new HMM has (r + 1)m states, each of which has the same degree as in M. We compute the total forward probabilities of all states (j, d) at position n, and this gives the probability of B(λ, r) for the generalized Hamming distance, in runtime r times that of the forward algorithm on the original sequence, or at most O(m²nr).

Hardness results

We might hope that we can compute the labelling for which B(λ, r) is of highest probability for a given r and one of our distance measures. Unfortunately, this is NP-hard, because the problem of computing the maximum probability labelling is also NP-hard [6], even for a fixed HMM with only two labels.

Optimizing the ball with radius 0, B(λ, 0), is exactly this problem, so optimizing balls for general radius r is clearly NP-hard for either the border shift or border shift sum distances. Similarly, it is NP-hard for the generalized Hamming distance, since if we used radius r = 0 and the label distance matrix H with value 1 for mismatched labels and 0 for matched labels, we are again maximizing the probability of a labelling. We conjecture that this problem stays NP-hard when restricted to the pure Hamming distance on state paths, and not labels.

Furthermore, maximizing the probability of a ball of labellings under any natural distance remains NP-hard even if we restrict attention to a single footprint f. The original HMM used in the proof that maximizing the label probability is NP-hard has only two labels, so there are only 2n possible footprints for a sequence of length n in this HMM. As such, if there existed a polynomial time algorithm that finds a maximum probability ball for a given footprint f, we could solve the most probable labelling problem by setting r = 0 and running that algorithm for each of the O(n) footprints.

Optimizing balls of states and paths: local search and special cases

Still, we can hope to find good balls of labellings using local search or global optimization in special cases. Here, we give a very efficient local search procedure to help optimize the probability of a ball of labellings using the border shift distance, and a global optimization procedure for computing the best centre for a ball of paths, not labellings, with the restriction that the path that forms the centre of the ball must divide the sequence up into long intervals, all using a single state.

Global optimization for paths with long intervals

While the general problem of finding the most probable ball of labellings is NP-hard, the problem can be solved in polynomial time for a special case. Here, we look for the most probable ball of state paths rather than labellings (which corresponds to each state having its own label). Furthermore, we are only interested in finding the most probable among those balls whose centre paths have each state lasting at least 2r + 1 positions. A path satisfying this requirement is called (2r + 1)-regular. We call a path π weakly (2r + 1)-regular if all its states except the final one last for at least 2r + 1 positions, while the final state lasts at least r + 1 positions. A weakly (2r + 1)-regular ball is a ball centered at such a path. Our dynamic programming algorithm maintains variables F [j, s] corresponding to the probability of the most probable weakly (2r + 1)-regular ball up to position j whose centre has state s at position j. The key idea behind the algorithm is that, by storing the probabilities of weakly (2r + 1)-regular balls, we can ensure that the last state of each path in the ball is the same. This allows us to devise a dynamic programming procedure similar to the Viterbi algorithm.

An optimal weakly (2r + 1)-regular ball up to position j can arise either by extending the final state of some centre path optimal to position j - 1 or by allowing a border to occur at position j - r because the last border must occur at least r positions previously. The recursive formula for calculating the probability of the new optimal ball is given by

where , the probability that the model emits the sequence y_j-2r-1to y _j in only the states s' and s, and with a forced transition between those states in the interval, conditioned on starting the interval in state s'; we write as a shorthand for x _i , x_i+1, ..., x _j = s, ..., s.

In this formula, the first term computes the probability of the best weakly (2r + 1)-regular ball ending at s at position j where there is a border at position j - r, and the second computes the probability of extending the previously optimal path in the same state.

Once we have computed the values of F [j, s] up to j = n - r, we multiply each of these values by a sequence of transition and emission probabilities corresponding to extending the state s up to position n. This ensures that our centres are (2r + 1)-regular in the strong sense. Finally, we pick the largest probability obtained in this way and reconstruct the actual state path.

The runtime of the algorithm is determined by the time needed to calculate the Q(s, s', j) values. We can calculate these sums in amortized constant time by reusing the values obtained for the previous value of j. This yields a runtime of O(nm²), as for the Viterbi algorithm.

Data used

We used the larger of the two data sets used to train Phobius [2]. The data set contains 247 transmembrane proteins, whose average length is 385 residues. The average number of transmembrane helices per protein is 5. We also used the Phobius HMM for our experiments. The training procedure for the model was originally designed for use with the 1-best algorithm, which is used in Phobius. For details about how the model was trained, see the original paper [2].

We used two error measures to assess the accuracy of our algorithms. The first counts a prediction as correct if its footprint is correct and the predicted location of each transmembrane segment overlaps with its true location in at least 5 positions. As we noted above, the authors of Phobius used this measure to identify correct predictions, so we call it the "Phobius" measure. The second measure we used counts a prediction as correct if its border shift distance to the true labelling is at most 5. We call this the "τ = 5" measure.

We note that the τ = 5 measure is much more stringent than the Phobius measure. In our experiments, we discovered far lower success rates for decoding algorithms according to this measure than to the Phobius measure. That measure does not penalize predictions that move the borders of a transmembrane segment far away from its centre and also allows more error when the prediction considers a border to be too close to the centre.

Initial labellings

We tested several different methods of choosing initial ball centres. For each of these methods, we then chose the ball of highest probability among those found by our local search procedure. In the first series of experiments, the initial ball centres were obtained by sampling from the conditional distribution of labellings given the protein sequence, using a standard algorithm [7, 8].

The second initialization method that we tested was choosing labellings corresponding to the K most probable state paths through the HMM. The K most probable paths for each protein were obtained using the software by Brown and Golod [9].

Finally, to see if we could use the prediction of Phobius as a good starting point for local optimization, we performed local search on the single labelling obtained from the 1-best algorithm in Phobius, and from K - 1 sampled paths augmented by the single labelling from Phobius.

Overall results

For this initial set of experiments, we used 10 samples from a variety of methods and optimized using our local search procedure.

Our methods are comparable to Phobius without local optimization for its measure of correctness (162/247 vs. 166/247). For our "τ = 5" measure, our optimized predictions from samples of the HMM space are significantly better (93/247 vs. 79/247). We used McNemar's test to verify that this result is significant at the 0.05 level for a 1-tail test.

In this experiment, our prototype implementation of our local optimization method required approximately 13 times as much time to compute predictions as did Phobius. However, we note that the parameters chosen (for r, for the number of iterations in the local search, and for the number of starting samples) were not done so as to reduce runtime, but to explore the effectiveness of our method.

Robustness of results

The behaviour of the local search algorithm depends on several parameters: the radius r of the ball whose probability we seek to maximize, the number of starting points, and the maximum number of local search iterations. We ran our sampling procedure for a variety of parameter choices. Our method is robust to these parameter choices: most experimental setups differed in accuracy by at most 2 - 3 proteins according to the Phobius measure. For the "τ = 5" measure, accuracies differed in a range of about 7 - 8 proteins. Even for the latter measure, the only parameter that noticeably affected the results was the ball radius r. Increasing the number of starting points from 10 to 20 had almost no influence on the Phobius measure accuracy of the local search procedure with sampling (results not shown). True footprints of high probability were very likely to be sampled at least once with only 10 samples. On the other hand, if the true footprint had low probability, the local search procedure would likely have chosen a ball of higher probability instead, and more samples would not help.

Changing the maximum number of local search iterations also did not influence the results significantly. We tested our algorithms with 6, 20 and 80 local search iterations; this had no effect on the Phobius measure, while letting the algorithm come closer to convergence improved the "τ = 5" measure at at most 2 proteins. Changing the radius of the ball had little influence on accuracy according to the Phobius measure (at most 2 for all tested values of r ≥ 2), but it did affect the performance according to the "τ = 5" measure. The number of correct predictions according to that measure varied from 84 (34%) to 93(37.6%) when we fixed the other parameters. Not surprisingly, the highest accuracy is achieved at r = 5. On the other hand, setting r = 6 yields the highest accuracy for the Phobius measure. Regardless of the radius used, the accuracies for the "τ = 5" measure were better than the accuracy obtained by the 1-best algorithm.

Experimental conditions

Unless otherwise noted, we used r = 5 in an attempt to directly optimize the "τ = 5" measure. We also show the results for r = 6 for sampling since this radius value yielded the best accuracy according to the Phobius measure. For each algorithm, we used 10 starting points for the local search, except the one where the local search was carried out from the 1-best prediction from Phobius.

Substituting one of the random samples with the 1-best prediction does not noticeably change the accuracy of the local search algorithm. Local searches starting from the 1-best prediction rarely produce higher probability balls than local searches based on sampling. Running local search from the Phobius prediction alone substantially increases the "τ = 5" accuracy for those predictions while not giving a noticeable performance gain according to the Phobius measure.

Initializing local search from labellings corresponding to the 10 most probable paths yields substantially worse accuracies than the other methods. This may be because of a low diversity of footprints among these paths.

Analysis of errors

The efficiency of our local search procedure when combined with sampling is not surprising. Since we are trying to optimize for high probability balls, the probability of sampling a labelling somewhere within the most probable ball should be large. As a result, the distance between a sample and the centre of the most probable ball will often be small, which makes the local search procedure likely to achieve at least a good approximation of the most probable ball. Obviously, this strategy will perform worse if the most probable ball has low probability, but in such cases we are not likely to obtain a good prediction anyway from any decoding strategy.