Sequence alignment by passing messages

A brief overview of pair hidden Markov models

The pair-HMM [1, 12] is a generative sequence model that can simultaneously generate a pair of aligned symbol sequences. This is different from the traditional HMMs, which generate only a single symbol sequence at a time [13]. Figure 1 shows two examples of pair-HMMs that are widely used in biological sequence analysis. As shown in Figure 1, a typical pair-HMM consists of three hidden states I_x, I_y, and M, which are used to model insertions in sequence x, insertions in sequence y, and matched (i.e., aligned) symbols in both sequences, respectively. The pair-HMM generates an aligned sequence pair (x, y) by making transitions between the hidden states according to the specified state transition probabilities. At state I_x, the model emits a symbol only to sequence x, while at I_y, a symbol is emitted only to sequence y. On the other hand, at state M, the model emits a pair of aligned symbols, where one symbol is added to x and the other symbol is added to y. Figure 1(C) gives an example of a sequence pair (x, y) that is generated by a pair-HMM. In this example, the underlying hidden state sequence that gives rise to the two sequences x = AACCG and y = CCGTT is I_xI_xMMMI_yI_y. This indicates that the first two symbols (i.e., AA) in x and the last two symbols in y (i.e., TT) are "insertions," which do not have any matching counterpart in the other sequence, while the last three symbols in x and the first three symbols in y (i.e., CCG in both sequences) are jointly generated by the pair-HMM, hence closely match each other. As we can see from this example, we can unambiguously identify the alignment of a given sequence pair (x, y), once the underlying hidden state sequence yielding the sequence pair is known. Of course, the hidden state sequence is generally not known, but there exist efficient algorithms that can be used for its prediction. For example, we can use the Viterbi algorithm [14] to predict the optimal hidden state sequence that maximizes the observation probability of the sequence pair (x, y). Alternatively, we can also predict the state sequence that maximizes the expected number of correct states, by first estimating the alignment probabilities between the symbols in x and y through the forward and backward procedures [13] and then applying the Needleman-Wunsch algorithm [15]. This will lead to the MEA alignment between the two sequences x and y.

Limitations of pair-HMMs

For the pair-HMM shown in Figure 1(A), any hidden state sequence $s=s_1{s}_2\cdots s_9$ such that s₁ = s₂ = s₃ = M and $s_4{s}_5\cdots s_9$ is a permutation of I_x I_x I_x I_y I_y I_y would lead to the sequence alignment shown in (1). When using this pair-HMM for predicting the optimal alignment of a sequence pair with the largest probability, this ambiguity may lead to performance degradation as these potential state sequences compete against each other. For this reason, it is generally more desirable to estimate the symbol alignment probabilities via the pair-HMM by considering all potential alignments and state sequences and use the estimated probabilities to find the MEA alignment that is expected to have the maximum number of correctly aligned symbols [9–11]. However, the aforementioned ambiguity also negatively affects the quality of the estimated symbol alignment probabilities, which is especially noticeable for sequence pairs with low percentage identity. In some cases, the alternative pair-HMM shown in Figure 1(B) is used to avoid such ambiguity. This alternative pair-HMM blocks transitions between the insertion states I_x and I_y, thereby prohibiting the model from inserting unaligned symbols to both sequences. For example, the alignment shown in (1) would not be allowed based on this alternative pair-HMM. However, due to this restriction, the pair-HMM in Figure 1(B) has a relatively stronger tendency to align unrelated sequence regions by treating them as mutations. This may again negatively affect the quality of the symbol alignment probabilities estimated based on the pair-HMM.

Another potential drawback of pair-HMMs is that the associated algorithms (i.e., the Viterbi, forward, and backward algorithms) can become numerically unstable for long sequences. Application of pair-HMMs to biological sequence analysis involves computing extremely small probabilities, which decrease exponentially with the sequence length. For example, based on the pair-HMM that was used in [9], the observation probability (i.e., the probability that the HMM may generate a given sequence pair) of a protein pair is typically in the order of 10^-230 for proteins of length 80, 10^-280 for proteins of length 100, and 10^-320 for proteins of length 120. As a result, pair-HMM algorithms are prone to underflow errors, unless they are carefully implemented to keep them numerically robust. So far, a number of schemes have been proposed to address this issue, such as using log transformations of the probabilities or normalizing the probabilities to keep them within a reasonable numerical range, and have been shown to work well for relatively long sequences [1]. However, log transformations can make the forward and backward algorithms considerably slower, and the normalization approach can still lead to underflow errors as the sequences get longer.

One further disadvantage of pair-HMMs is that the algorithms that are used with the model cannot be easily parallelized. Although the Viterbi, forward, and backward algorithms for pair-HMMs are relatively efficient, they are still computationally expensive to be used with very long sequences. Moreover, as the algorithms are not amenable to massive parallelization, this makes the pair-HMMs not suitable for large-scale sequence analysis tasks, such as the whole genome alignment, despite their superior performance compared to other heuristic methods.

A message-passing scheme for estimating symbol alignment confidence scores

Here, we propose a novel method for aligning biological sequences that can effectively address the aforementioned shortcomings of pair-HMMs. The proposed method is based on a message-passing scheme, where messages are iteratively exchanged between neighboring symbol pairs to estimate the level of confidence for potential pairwise symbol alignments. The main underlying motivation is to develop an "analytical" method that can directly estimate the symbol alignment probabilities, without specifically modeling symbol insertions and deletions. This stands in contrast to the pair-HMM approach, which is essentially based on a "generative" sequence model that tries to explicitly model symbol insertions/deletions, in addition to symbol alignments. As discussed before, modeling symbol insertions in pair-HMMs can lead to subtle issues with potentially negative effects, and considering that our ultimate goal lies in finding an accurate sequence alignment through effective estimation of the symbol alignment probabilities, a method that can directly estimate these probabilities without explicitly modeling insertions/deletions would be desirable.

Suppose x = x₁x₂ ⋯ x _L and y = y₁y₂ ⋯ y _M are the two sequences to be aligned. We define c _xy (i, j) as the symbol alignment confidence score between x _i (the i-th symbol in x) and y _j (the j-th symbol in y). The score c _xy (i, j) provides a quantitative measure of confidence as to whether x _i and y _j should be aligned to each other or not, and we assume c _xy (i, j) ∝ P(x _i ~ y _j |x, y), where P (x _i ~ y _j |x, y) is the posterior symbol alignment probability between x _i and y _j given the sequences x and y. We estimate the alignment confidence score by iteratively passing messages between neighboring symbol pairs, where each symbol pair (x _i , y _j ) corresponds to a potential symbol alignment in the true (unknown) sequence alignment between x and y. For example, during the estimation process, the symbol pair (x _i , y _j ) will exchange messages with its two neighbors (x_i-1, y_j-1) and (x_i+1, y_j+1), and similarly, the pair (x_i+1, y_j+1) will exchange messages with (x _i , y _j ) and (x_i+2, y_j+2). The message-passing process is illustrated in Figure 2, where the solid lines indicate the messages that are used to update the alignment confidence score c _xy (i, j) of the symbol pair (x _i , y _j ). The dashed lines correspond to messages that are used to update the confidence scores of other symbol pairs.

The pseudocode of the proposed message-passing algorithm is as follows:

   STEP-1 Initialize c _xy (i, j).

   STEP-3 Normalize c _xy (i, j).

   STEP-4 If c _xy (i, j) has converged, then terminate the algorithm.

                  Otherwise, go to STEP-2.

In STEP-1, we first initialize the alignment confidence score c _xy (i, j), where we can simply use random initialization. If a preliminary sequence alignment of x and y is available (e.g., obtained from a simple heuristic method), we can also initialize the score based on this alignment such that c _xy (i, j) = 1 if x _i and y _j are aligned, and c _xy (i, j) = 0 otherwise. Next, in STEP-2, the alignment confidence score c _xy (i, j) of the symbol pair (x _i , y _j ) is updated based on the scores of its two neighbors (x_i-1, y_j-1) and (x_i+1, y_j+1). Note that the score is set to c _xy (i, j) = 0 if i ∉ {1, ⋯ , L} or j ∉ {1, ⋯ , M}. P (x _i , y _j ) is the joint occurrence probability of the symbol pair (x _i , y _j ), which is essentially equivalent to the joint emission probability of an aligned symbol pair (x _i , y _j ) at the match state M of a pair-HMM. It should be noted that this probability P (x _i , y _j ) is not location-dependent and is simply determined by the symbols x _i and y _j . The weight parameter λ ∈ [0, 1] is used to balance the contribution from the neighbors and that from the joint probability of (x _i , y _j ) in estimating the alignment confidence score. A large λ gives more weight to the "messages" received from the neighbors in estimating the scores, which tends to penalize gaps more heavily, and it generally leads to longer aligned regions with fewer gaps. On the contrary, a small λ gives more weight to the joint symbol occurrence probability P (x _i , y _j ) while giving less weight to the messages received from the neighbors, which tends to be more lenient to gaps. Once the symbol alignment confidence score c _xy (i, j) is updated for all i = 1, ⋯ , L and j = 1, ⋯ , M, we normalize the scores to keep them within a proper numerical range, as shown in STEP-3. For example, a simple way would be to divide the score matrix C = [(c _xy (i, j)] by its matrix norm to normalize the confidence scores. After normalization, the updated scores are compared to the scores in the last iteration, and the algorithm terminates if the specified convergence criterion has been met. Otherwise, the algorithm goes back to STEP-2 and repeats the message-passing process.

Dataset and experimental set-up

In order to evaluate the performance of the proposed message-passing scheme, we carried out pairwise sequence alignment experiments based on the BAliBASE 3.0 protein alignment benchmark [16]. BAliBASE is arguably the most widely used benchmark for multiple sequence alignment, and it has been utilized by most multiple sequence alignment algorithms for assessing their performance. The benchmark consists of five reference sets, where Reference 1 consists of two subsets: V1 and V2. Each reference set consists of multiple sequence alignments that satisfy specific criteria, such that different reference sets can be used to test the performance of multiple sequence alignment algorithms under different conditions. For example, each alignment in Reference 2 consists of sequences that share reasonably high identity (> 40%) and "orphan sequences" that share little identity (< 20%) to other sequences in the alignment. Reference sets 4 and 5 are constructed such that every sequence has at least one other sequence in the same alignment whose identity exceeds 20%. Sequences in Reference 4 and Reference 5 may contain large N/C-terminal extensions or internal insertions, respectively. Further details of the BAliBASE 3.0 benchmark can be found in [16].

For every sequence family in BAliBASE 3.0, we performed pairwise sequence alignment for all possible sequence pairs in the given family. The pairwise alignment was performed in the following manner. First, we estimated the probabilistic symbol alignment confidence score using the proposed message-passing scheme. In our experiments, we used three different values of λ (= 0.25, 0.5, and 0.75) to investigate the effect of λ on the overall sequence alignment performance. For the joint symbol occurrence probability P (x _i , y _j ), we used the joint emission probability (at state M) of the pair-HMM that was used in [9]. At the end of each iteration, we normalized the alignment confidence score by dividing the confidence score matrix C by the matrix 2-norm: C ← C /|| C ||₂. We terminated the message-passing process if ${\sum }_i{\sum }_j|c_{xy}\left(i,j\right)-{\tilde{c}}_{xy}\left(i,j\right)|<0.01$, where c _xy (i, j) is the current score and ${\tilde{c}}_{xy}\left(i,j\right)$ is the score obtained in the previous iteration. Once the scores converged, based on our assumption that c _xy (i, j) ∝ P (x _i ~ y _j | x, y), we used the confidence score c _xy (i, j) to find the MEA alignment through dynamic programming. The predicted alignment was compared to the benchmark alignment in BAliBASE 3.0 to compute the sensitivity $\left(\mathsf{\text{SN}}\right)=\frac{TP}{TP+FN}$ and the positive predictive value $\left(\mathsf{\text{PPV}}\right)=\frac{TP}{TP+FP}$, where TP is the number of correctly aligned symbol pairs, FP is the number of incorrectly aligned pairs, and FN is the number of symbol pairs that are aligned in the benchmark alignment but not aligned in the predicted alignment. For comparison, we repeated similar experiments using the pair-HMM with the same set of parameters as the one used in [9].

Performance of the proposed message-passing scheme

From Table 1, we can clearly see that the proposed message-passing scheme significantly outperforms the pair-HMM approach in terms of SN and PPV, for all three values of λ. For example, the message-passing scheme achieved up to 0.23 higher SN and 0.09 higher PPV for λ = 0.25, and up to 0.37 higher SN and 0.19 higher PPV for λ = 0.75. Our experiments showed that a larger λ tends to yield more accurate alignments, while a smaller λ tends to make the algorithm converge faster, hence computationally more efficient. For example, when the weight parameter was set to λ = 0.25, the message-passing scheme was around 2.3 ~ 2.5 times faster than the pair-HMM, while still yielding much more accurate alignments.

The results in Table 1 demonstrate that, on average, the proposed message-passing scheme considerably improves the quality of sequence alignment over the traditional pair-HMM approach. In order to see whether the proposed scheme also leads to a consistent improvement for most sequence pairs, we calculated the difference between SN _MP (the sensitivity of the message-passing scheme) and SN _HMM (the sensitivity of the pair-HMM-based approach) for every pairwise sequence alignment that we have performed in our experiments. Similarly, we calculated the difference between PPV _MP (the PPV of the message-passing scheme) and PPV _HMM (the sensitivity of the pair-HMM approach) for all sequence pairs in BAliBASE 3.0. Figure 3 shows the distributions of SN _MP - SN _HMM and PPV _MP - PPV _HMM for all sequence pairs. To avoid any bias from unsuccessful alignments, sequence pairs for which neither method yielded an alignment with at least one correct symbol alignment were excluded. The plots in the left column of Figure 3 show the distributions of SN _MP - SN _HMM , and those in the right column show the distributions of PPV _MP - PPV _HMM . The results obtained from the same reference set are shown in the same row, where the first row shows the results on RV11 and the last row shows the results on RV50. As we can see in Figure 3, every single distribution shown in the figure has a much larger probability mass in the right-half plane, which clearly demonstrates that the proposed message-passing scheme consistently outperforms the pair-HMM-based approach for most (though not all) sequence pairs. In many cases, the improvements in SN and PPV were quite significant (0.4 ~ 0.8), which shows that the proposed scheme can often find an accurate sequence alignment even when the pair-HMM has difficulty aligning the sequences.