Hidden Markov models of biological primary sequence information ( multiple sequence algnments / protein modeling / adaptive algorithms / sequence classification )

Hidden Markov model (HMM) techniques are used to model families of biological sequences. A smooth and convergent algorithm is introduced to iteratively adapt the transition and emission parameters of the models from the examples in a given family. The HMM approach is applied to three protein families: globins, immunoglobulins, and kinases. In all cases, the models derived capture the important statistical characteristics of the family and can be used for a number of tasks, including multiple alignments, motif detection, and classification. For K sequences of average length N, this approach yields an effective multiple-alignment algorithm which requires O(KN2) operations, linear in the number of sequences. Comparative analysis of primary sequence information is a major tool in the elucidation of the molecular mechanisms of replication and evolution of organisms and the structure and function of proteins. For the simple case of pairwise sequence comparison, good algorithms exist (see refs. 1 and 2 for recent reviews) that can align two sequences of length N in roughly O(N2) steps. Most of these algorithms are based on dynamic programming (3), with location-independent substitution and gap penalties. Unfortunately, when dynamic programming is applied to a family of K sequences its behavior scales like O(NK), exponentially in the number of sequences (4). A number of algorithms have been devised to try to tackle the multiple alignment problem (see refs. 5-7 for some of the most recent ones). Most protein sequence relationships exhibiting >50%o identical residues can be aligned by several of these algorithms. Many of the most interesting protein families, however, exhibit conservation far below 50%o identity. To date, alignment methods have not been developed that can correctly identify all the motifs that define each protein family (2). Here, we apply a different approach, based on hidden Markov models (HMMs), to the problem of modeling and aligning a family by using primary structure information only. Initial results were presented (8). Markov models and the related expectation-maximization (EM) (9) algorithm in statistics have already been applied to biocomputational problems (10-13). Krogh et al. (14) were the first to demonstrate the power of a similar method on the globin family. Rather than starting from pairwise alignments, the approach seeks to take advantage of the massive amount of information typically present in a family with a flexible use of positiondependent parameters. A new algorithm is introduced for the iterative adjustments of the parameters of the models. The algorithm is used here to model three protein families: globins, immunoglobulins, and kinases.tt HMMs and Learning A first-order discrete HMM (15) is completely defined by a set of states S, an alphabet of m symbols, a probability transition matrix T = (tv), and a probability emission matrix E = (eta). When the system is in state i, it has a probability t(/ of moving to state] and a probability eia of emitting symbol a. Only the output string is observed, one of the goals being the reconstruction of the underlying hidden transitions. As in the application of HMMs to speech recognition, a family ofbiological sequences can be seen as a set ofdifferent utterances of the same word generated by a common underlying HMM with a left-right architecture, with m = 4 for DNA or RNA and m = 20 for proteins. Common knowledge about evolutionary mechanisms suggests to introduce three classes of states (in addition to the start and end states): the main states, the delete states, and the insert states with S = start, ml,. . . ,MN, il,. . . ,JN+l, di,. . . ,dN+l, end . N is the length of the model. The main and insert states always emit a letter of the alphabet, whereas the delete states are mute. The linear sequence of main state transitions is the backbone of the model. Self loops on the insert states allow for multiple insertions. Architectural variations are possible and may be tailored to particular problems when additional information is available. Given a set of training sequences, the parameters of a model can be iteratively modified to optimize the fit of the model to the data according to some measure, usually the product of the likelihoods of the sequences. Different algorithms are available forHMM training, including the classical Baum-Welch algorithm (15). Here, we introduce a smooth algorithm which is particularly simple and can be used on-line-i.e., after the presentation of each example. The mathematical properties of this algorithm and its relation to other approaches have been studied (16). We first reparametrize the model by using a new set of variables wV and via in the form tu, = eWu/yk e wi and ej, = evi/: eve This reparametrization has two advantages: (a) modification of the w and v parameters automatically preserves the normalization constraints on probability distributions and (b) transition and emission probabilities can never reach the absorbing value 0. We then iteratively cycle through the training set and compute, for each sequence, the corresponding most likely path through the model. This can be done efficiently in O(N2) steps by using a dynamic programming scheme known as the Viterbi algorithm. Being in a state i along a Viterbi path, we update the parameters of the model according to AwV = q(Tv tY) and Avia = 'q(Eia eia), where v1 is the learning rate.  HMM, hidden Markov model. Present address: Department of Biological Sciences, University of Nevada, Las Vegas, NV 89154. ItAIl programs, technical reports, data sets, complete alignments, and other results are available, for noncommercial purposes, from the authors upon request. 1059 The publication costs of this article were defrayed in part by page charge payment. This article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. §1734 solely to indicate this fact.  Side states d, (resp. i,) correspond to deletions (resp. insertions). each step of a Viterbi path, and for any state i on the path, TV = 1 (resp. Eta = 1) if the i --j transition (resp. emission of a from i) is used and 0 otherwise. Parameters of a loop are updated every time the loop is traversed. These rules are repeated for each example until no significant variations occur. This algorithm approximates a gradient descent procedure on the negative log-likelihood of the sequences (16). As such, it can be expected to converge to a possibly local maximum likelihood estimator. Once aHMM has been successfully trained, it can be used for a number oftasks. Multiple alignments result immediately from aligning the corresponding optimal paths. Discrimination ofwhether any given sequence belongs to the family can Proc. Nati.  (1994) be based on its likelihood according to the model. Structural properties ofthe model revealed-for instance, by an entropy plot ofthe emission distribution along the backbone-are also useful.  Experiments have been performed with several protein families, including globins, immunoglobulins, kinases, G-protein-coupled receptors, aspartic proteases, and human immunodeficiency virus membrane proteins. The focus here will be on globins, immunoglobulins, and kinases. In all three experiments, the length of the model is initially set to the average length of the sequences in the family and then kept constant for the entire experiment. The probability parameters of the models are initialized uniformly prior to training. A random subset of sequences is typically used for training with a learning rate v1 = 0.1. The remaiing sequences are used for validation. After training, all the sequences in the training and validation set are aligned to the model. For brevity, only aportion of the resulting multiple alignments is displayed, by using a subset of phylogenetically representative sequences covering both the training and validation sets. Globins. The globins form a well-known family of hemecontaining proteins that reversibly bind oxygen and are involved in its storage and transport. From crystallographic studies, all globins have similar three-dimensional structure characterized by seven a-helices, labeled A, B, C, E, F', G, and H (some structures also have a short D helix) (17). The globin sequences used here were extracted from the nonredundant (NR) data base composed of Protein Identification Resource 34.0, Swiss-Prot 23, and GenPept (translated Gen-