Hidden Markov Models in Computational Biology

zyxwvutsrqpo zyxwvutsrq zyxwvutsrq J . Mol. Bioi. (1994) 235, 1501-1531 zyxw zyxwvu Hidden Markov Models in Computational Biology Applications to Protein Modeling Anders Krogh'f, Michael Brown', I. Saira Mian' Kiminen Sjolander' and David Hausder'S zyxwv 'Computer and Information Sciences 2Sinsheimer Laboratories University of California, Santa Cruz, C A 95064, U.S.A. . Ridden.Markov Models (HMMs)are applied t.0 the problems of statistical modeling, database searching andmultiple sequencealignment of protein families andprotein domains. These methods are demonstrated the on globin family, the protein kinase catalytic domain, and the EF-hand calcium binding motif. In each case the parameters of an HMM are estimated from a training set of unaligned sequences. Afterthe HMM is built, itis used to obtain a multiple alignment of all the training sequences. It is also used to search the SWISS-PROT 22 databaseforother sequences. thatare members of the givenprotein family, or contain the givendomain. The Hi" producesmultiplealignments of good quality that agree closely with the alignmentsproduced by programsthat incorporate threedimensional structural information. When employed in discrimination tests (by examining howclosely the sequences in a database fit the globin, kinase andEF-hand HMMs), the HMM is able todistinguish members of these families from non-members with a high degree of accuracy.Both the HMM and PROFILESEARCH ( a technique used t o searchfor relationships between a protein sequence and multiply aligned sequences) perform better in these tests than PROSITE (a dictionary of sites and patterns in proteins). The HMM appecvs t o have a slight advantageover PROFILESEARCH in terms of lower rates of false negativesand false positives, even thoughthe HMM is trainedusing onlyunaligned sequences, whereas PROFILESEARCH requires aligned training sequences. Our results suggest the presence of an EF-hand calciumbindingmotif in a highlyconserved and evolutionary preserved putative intracellular region of 155 residues in the a-l subunit of L-type calcium channels which play an important role in excitation-contraction coupling. This region has beensuggested to contain the functional domains that are typical or essential for all L-type calcium channels regardless of whether they couple to ryanodine receptors, conduct ions or both. zyxwv '\ \ ' Keywords: hidden Markov models; multiple sequence alignments; globin; kinase; EF-hand 1. Introduction The rate of generation of sequence data in recent years providesabundantopportunitiesforthe development of new approachesto problems in computational biology. Inthispaper, we apply - t Present address: Electronics Institute, Build 349, Technical University of Denmark, 2800 Lyngby. Denmark. $ Author to whom all correspondence should be oddreased. 8 Abbreviations ' HMM, hidden Marko,. ,,,dels; m,Expectetion-Maximization;ML. maximum like€ihood;MAP, maximum a posten'wi; hTL-score. negative log likelihood score. CO22-2836/94!031301-31 $08.00!0 hidden Markov models (HMMdl to the DrobiemS of statistical modeling, 'database searihing, and multiple alignment of protein families and protein domains. To demonstrate the method, we examine three protein families. Each family consists of a set of proteins that have the sameoveraIl three-dimensional structurebut widelydivergent. sequences. Features of the sequences that are determinants of folding, structure and functionshould be present as conserved elements in the family of sequences. We consider t.he ranging globins, proteinswhole in lengthfrom 130 t o 170 residues (with few exceptions) and two domains, the protein kinase catalytic domain (250 to 300 residues)and theEF-hand caicium-binding motif (29 residues). The same 1.501 'r; 199-1.-\cdcrnir I'rpss Limited 1502 zyxwvutsrqpo zyxwv zyxwvut zyxwvutsrq zyxwvuts IIidden JZarkoz: :Models approach can beused to model families of nucleic acid sequences as well (Krogh et al.. 19936). A hidden Markor model (Rabiner, 1989) describes a series of observations by a "hidden'' stochastic process, a Markov process. In speech recognition. where HMMs hare been used extensively, the observationsaresounds forming a word. anda model is one that by i t s hidden random process generates these sounds with high probability. Every possible sound sequence can be generated b.v the model with some probability. Thus: the model defines a probabiIity distribution over possible sound sequences. X good word model would assign high probability t o all sound sequences t.hatare likely utterances of the word it models. and low probability to any othersequence. In this paperwe propose a n HM31 similar to the ones used in speech recognition t o model protein families such as globins and kinases. I n speech recognition. the '-alphabet" from which words are constructed could be the set of phonemes valid foraparticular language: in protein modeling. thealphabet we use is the 20 aminoacids from which protein molecules are constructed. t7'here t.he observations in speech recognit.ion are words. or strings of phonemes. in protein modeling the observations arestrings of a.mino acids forming the primarysequence of a protein. A model for a set of proteins is one that assignshighprobabilitytothe sequences in that particular set. The HYY w e build identifies a set of positions 'fhat describe the (more or less) conserved first-order &ructure in the sequences from a given family of proteins.In biological terms.thiscorresponds to identifl-ing the core elements of homologous molecules. The model provides additional information, such as the probability of initiating an insertion at any position in the model and the probability of extendingit. The structure of the madel is similar to that of a profile (Waterman & Perlu-itz, 1986: Barton Q Sternberg, 1990: (;rihsko\et al., 1990: Bowie et ai.. 1991: Iiithy et nl.. 1991). but slightlv more general. Onw we have built the model from unaligned sequences. we can generate a multiple alignment of the sequences using a dynamic programming method. By employing it for database searching, the modelcan be used to tliscriminate sequences that belong to a given family from non-members. Finally. we can study themodel we havefound directly. and see what it reveals about the common structure underlying the various sequences in the famill-. Our method of multiple alignment differs quite markedly from conventional techniques. which are usually based on painvise alignments generated by d y n m i c programming. schemes (Waterman. 1989: Fend & DootittIe, 1987; Barton, 1990: Subbiah & Harrison, 1989). The alignments produced by these methodsoften depend strongl? on theparticular vaiues of t h e pa.rameters required by the method. in particular the gap penalties (Vingron & Argos. 1991). Furthermore, a given set of sequences is likely to possess both fairly conserved regions and highly variable regions, vet conventional global methods assign identical penalties for all of the sequences. Substitutions, insertions,or deletions in a regionofhigh conservation should ideally be penalized more than in a variable region, and Some kinds of substitutionsshouldbe penalized differently in one position compared t o another. That is one of the motivations for the present work. The statistical model we propose corresponds to multiple alignment with variable. position-dependent gap penalties. Furthermore, these penalties are in large part learned from the data itself. Essentially, we build astatist.ical model duringthe process of multiple alignment. rather than leal-ing this as a separat.e task to be doneafter the alignment is completed. \Ye believe the model should guide the alignment as much as the alignment determines the model. We are not the first group to employ hidden Markov models i n computationa.i biology. Lander & Green (198T) used hidden Markov models in the construction of genetic linkage maps. Ot.her work employed HJIUs to distinguishcoding from noncoding regions in DSA (Churchill, 1989). Later, simple HJIJIs were used in conjunction with the EM algorithnl to model certain protein-binding sites in DSA (Lawrence Q R.eilly, 1990; Cardon & Stormo, 1992) and. more recently, to model the S-caps and C'-caps of alpha helices inproteins (D.Morris, unpublished results). These applications of HMNs and the EJl (Espertation-Maximization) algorithml including our own. presage a more widespread use of this technique in rotnputational biology. During the time that we have been developing this approach. several related efforts hare come to our attention. One is that of \illite. Stu1t.z and Smith (Whiteet nl.. 1991: Stnltz et d.. 1993), who use HMMs to model protein superfamilies. This work is more ambitious to than our o w l . since superfamiliesareharder characterize than families. I t is not get dear how successful their work has been since no results are reported for sequences not in the trainingset. If there are weaknesses in their method. it is possible that these are due to the use of handcrafted models and wliance on prealigned data for parameter estimation. T n contrast. our models have a simple regular structure. and we are ablet o estimate all the parameters of these models. including the size of the model directly from unaligned training sequences. Interestinglyenough,theyindependently propose an alternate H l i N state structure similar to oursi i n section 6.3 of theirpaper (White et al.. 1991). where they discuss the relationship of their work to Bowie and co-workers (Bowie et al., 1991): but they do not pursue this further. Tt is possible t.hat the txpe of models we use may work better for characterizing superfamilies thanthose investigated by \Vhite et 01. However. it is more likely that they are too simple. and that richer and more varied state zyxwvuts zyxwvu Hidden Markov Models zyxwvu zy 1503 zyxwvut zyxwvu zyxwvuts zyxwvuts Figure 1. The model. structure along the lines they propose is required for this problem. We recentlyfound that Asai et al. (1993) ha.ve applied HMMs to theproblem of predictingthesecondarystructure of proteins, obtaining prediction rates that are competitive with previous methods in some cases. In addition, Tanaka et al. (1993) also discuss the relationship between t h e HMM method for obtaining multiple alignments and previous methods. Finally. in work most closely related to our own, since the time ne presented a preliminary report on this work (Haussler Br Krogh, 1992; see also Haossler el ai.. 1992): BaIdi el ai. (1993) have further demonstrated the usefulness of this technique by producing multiple aiignments for immunoglobins and protease as we11 as globins and kinasest. 2. Methods (a) H Y , W architecture - zyxwvuts zyxwvu zyxw Consider a family of protein sequences that all have a commonthree-dimensionalstructure.forrsample the globins. The common structurein these sequences can be defined as a sequence of positions in spacu where amino acidsoccur. In the case of globins.whosr structure contains principally z-helices. the 1% or so helical positions have been named A I . -12. . ... A 16. HI. ... etc.. where the letterdenotesthen-helis.andthenumber indicates the location within that z-helix (see for esample Bashford et 01.. 1985). For each of these positions there is a (distinct)probabilitydistribution over the 20 amino acids that measures the likelihood of each aminoacid occurring in that position in a typiral g1,obin. as well as the p m b a b i l e that there is no amino acid in that position (Le. that a sequence belonging to this family may haveagapatthat position in a multiplealignment). These have been called profiles (Waterman 8: Perlwitz. 1986; Barton 8. Sternberg. 1990: Gribskov cl 01.. 19% Bowie et el.. 1991; Liithy et al.: 1991). A profile of globins can be thought of as a statistical model for the family of globins, in that for any sequenceof amino aeids. it defines a probability for that sequence:in such a way that globin sequences tend to have much higher probabilities than non-globin sequences. The type of hidden Markov model we use as a statistical model for a protein family can be viewed as a generalized profile. However. instead of describing the HM31 directly . .. t They have dere1op.I a variant i n terms of theprobability it assignstoeach protein sequence. we find that it is easier to first think of an HMM asastructurethatgenerates proteinsequences bya random process. This structure and corresponding random process is illustrated in Figure 1 and can be described as follo~vs. The main line of the HMM contains a sequence of Y states. \vhich we call match states, corresponding t o psitions in a protein or columns in a multiple alignment (M equals 1 in Fig. I ). Each of these states can generate a . letter x from the %letter amino acid alphabet. according 9(.rlmk). k = 1 ... N . Thenotation toadistribution 9 (.rlmt) means that each of the match states m,. 1 S k S .II.h a w distinct distributions. For each match state mk, there is a delete state dk thatdoes not produce anv amino acid but is a '-dummy" state used to skip mk. Finally, there BIP a total of .I1 + 1 insert states to either sideof tbe match states which generate amino acids in exactly tb same way as the match states. but use probability distributions 9 ( . r l i k ) . In Figure 1. match.deleteandinsert states are shown as boxes. circles and diamonds, respectively. For convenience, we have added a dummy -'HE(;IS" state and a dummy "EXD" state, denoted mo and mM-I . respectively. which do not produce any amino acid. From each state. there are three possible transitions to otherstates.alsoshown in Figure 1. Transitionsinto match or delete states always move forward in the model, whereas transitions into insert states do not. Sote that nlultiple insertions between match states can occur, since the self-loop on the insert state allows a transition from the insert state to itself. The transition probability from state q to stater is called 9(rlq).Our notation is summarized in Table 1. A sequencecan begenerated by a "random walk" through the model as follows: Commencing at state mo (BEGIS). choose a transition to m,. d,. or io randomly of the n~rthod described here that employs a gradient cirzwnt training algorithm in plare of the EM algorithm. Table 1 Notation X 3 L 9. tmlh .Y .I1 m. i. d ]no. n l y + , B(.elg) g(rlql Amino acid Sequence of amino acids (s= rl ...xL) Length of sequence Stnte in HMM A sequence of states. q, ...qN Somber of states in a yat.h Length of model Jktrh. insert, and delete states Begin and end s t a t e Probability distribution of amino acids in state q Probability of a transition from sbatr q to r 1504 zyxwvutsrqpo zyxwvu zyxwvut zyxwvutsr zyxwvut zyxwvut zyxwvutsrqp zyxwvutsr zyxwvutsrqpon Hidden Markov Models according to the probabilities 9 ( m , l m o ) , 9(dllmo), and 9(iolmo). If m, is chosen, generate thefirst amino acid rl from the probability distribution b(z)m,), andchoose a transition to the next state according to probabilities F (-]ml), where* indicates anypossible next state. If this next state is the insert state i,, then generate aminoacid zzfrom B(zli,) and select the next state from 3-(-lil). If delete Id2) is chosen next, generate no amino acid, and choose thenext state from Y(-ld2). Continueinthis to the END state, generatinga manneralltheway sequence of amino acids zl. z2... z, by following a path of states qo. ql . .. qN, qN+l through the model, where go = m, (the BEGIN state) and qN+I = mM+l (the Eh’D state).Becausethedelete states donot produceany amino mid,N is larger than or equal to L. If qi is a match or i n s e r t state, we define I(i) to be the index in the sequence x l ... zLof the aminoacid produced in state pi. The probability of the event thatt.he path qo .., q N + ,is taken and the sequencerl ...x, is generated is principle find the model t h a t best describes a given set of sequences. Given a set of training sequencesS( l ) , ...,s ( n ) ,one can see how well a model fits them by calculating the probability thatit generates them. Thisprobabi1it.y is simply a product of terms of the form given by equation(2), i.e. n Prob(sU)fmodel), n Prob(sequenceslmode1) = j- 1 (3) .. where each term Prob(s(j)l model) is calcul&ed by substituting z1 .. . xL = ~ ( j in ) equation (2). This is called the likelihood of the model. One would‘like this value to be high. The maximumlikelihood (PIIL) methodofmodel estimation is to find the mode1 that maximizes the likelihood (3). An alternate approach to ML estimation is the maximum’ a posferiuri (MAP) approach.Here, we assume a prior probability d h i b u t i o n over all possible parameters of the model embodying prior beliefs on what a model Prob(z, ... q , qo ... qn+ I model) should be like. This can thenbe used to “penalize” models that are known to be bad or uninteresting. We discussthis N MAP estimation?we t r y further in Rrogh et ul. (1993~). In =F(mN+IIqN) X n~(qiIqi-l)B(~,~i~Iqi), (1) i- 1 to maximize the posterior probability of the mode1 given the sequences. Using Bayes rule. the posterior probability where we set B(zIMlqi) = 1 if qi is a delete state. The can be calculated as probabiiity of a n y sequence zl ... tLof amino acids is a sumoverall possible pathsthat couldproduce that Prob(mode1lsequences) sequence. which we write as follows: Prob(mode1) - Prob(sequenceslmode1) Prob(z, ... zLlmodel) Prob(sequences) . (4) = Prob(z, ... z,, qo ,..‘qN+llmcdel). (2) Here Prob(mode1) is the prior probability distribution, P . l h r l 0 ...l N S 1 and Prob(sequences) can be viewed as a normalizing constant. Since this normalizing constant is independent In this way a probability distribution on the space of of the model. MAP estimation is equivalent to pquences is defined. The goal is to find a model (Le. a maximizing &roper model length andprobabilityparameters)that accuratelydescribes a Eamilyof proteinsby assigning Prob(sequencelmode1) Prob(mode1). t5) large probabilities to sequences in that family. Thisparticularstructureforthe HMM was chosen over all possiblemodels. The MAP approach isclosely because i t is the simplest model that captures the strucrelated to minimum description length (Jurka & tural intuition of a protein: (a) a sequence of positions. Milosa\-ljevic. 1991) and minimum mesage length each with i t s own distribution over the amino acids: (b) .(Allison el al.. 1992) methods. the possibility for either skipping a position or inserting There is no known efficientway to directly calculate the extra amino acids between consecutive positions; and (c) YL or MAP sense, best HMM model eitherinthe allowing For the possibility that continuing an insertion or However. there are algorithms that given an arbitrary deletion is morelikely thanstarting one. This choice startingpoint find a local maximum by iteratively appears to haveworked well for modeling t.he protein re-estimating the model in such a way that thelikelihood famiiies that we hsve examined, but othertypes of HMMs (or the posterior probabi1it.y) increases in each iteration. may be better at other tasks (e.g. the more elaborate The mostcommon one is the Baum-Welch or forwardmodels for proteinsuperfamilies used by White et ai.. backward algorithm (Rabiner. 1989; Lawrence & Reilly, 1991; Stultz el al., 1993). The important feature of the 1990), which is a version of the general EM method often HMM method is i t s generality. One can chooseany strucused in statistics (Dempster et al.. 1977). The process of ture for the states and transitions that is appropriate for the EM algorithm can be viewed as an iterative adapthe problem at hand. Examples of more general HMIY tation of the model to fit the trainingsequences. The steps architectures are given in sections (d) and (e). below. in this process can be summarized as follows: zyxwvutsrq zyxwvuts (b) Estimating the parameters of an H M N f r o m training sequences AI1 theparameters in the HMM (i.e. the transition probabilities and the amino acid distributions) could in principle be chosen by hand. from an existing alignmentof protein sequences,as in Gribskov et al. (1990), White el al. (1991).Stultz et al. (1993), or from information about the threedimensionalstructure of proteins, as inBowie et al. (19911, White et al. (1991), Stultz et aZ. (1993). The novel approachwetake is to “learn” the parameters entirely automatically from a set of unalignedprimary sequences, using a n EM algorithm. This approach can in (1) An initial model is created by assigning values to the transitionprobability Y (rip) andtheaminoacid generation probability P (xlq)for each2, q and r, where z is one of the 20 amino acids andq and r are states in the HMN connected by a transition arc. If one already knows some features present in the sequences, or constraints on the sequences,these maysometimes be encoded in the initial model. Thecurrent model is set to thisinitial model. (2) Csing the current model, all possible paths for each training sequence are considered in order to get a new estimate J (514) of-the transition probability 9( r l g ) and a new estimate b(zlq) of theamino acidgeneration probability P ( z ( q ) for each z. q and r. The transition zyxwvu zyx zy zyxwvuts zyxwvuts Hidden Markov Models zyxwvuts zyxwvutsrqp zyxwvutsr probability estimate g ( r l q ) is obtained by counting the number of times a transitionis made from state q to 7, for all paths of all training sequences, wsighted by the probP(zIq) is made in a ability of the path. The estimate of times the similarmanner.bycountingthenumber amino acid x is aligned to the stateq. (3) In the next step of ML estimation, a new current model is created- by simply replacing Y (714) by Y ( r l q ) and P(z1q) by b(rlq)for e3ch 2, q an! r. In MAP EM estimation, the parameters9 (slq) and d ( z J qare ) further modified by considering the prior probabilityof the model before they are used to replace the old parameters. ( 4 ) Steps (2) and (3) are repeated until t.he parameters of the current model change only insignificantly. . Since the quality of the current model (as measured by equations (3} or (5)) incre+ses in each iteration, and no model is arbitrarily good, the process eventually terminates and produces a model that is, at least locally, the best model for thetrainingsequences t.o withinsome specified precision of theparameters(Dempster el a.1.. 1975). Typically. this occurs very rapidly(e.g. in less than 10 iterations) even for large modelsand large sets of training sequences. The main computzitional bottleneck in the algorithm is step (2). since individualiy examining each possible path for every training sequence would generally take time exponential in the length of the longest training sequence. .However, it is possible to useadynamicprogramming technique known as the forward-backward procedure to speed up thisstep. Using this method, thenew parameter t,he estimates can be calculated in time proportional to number of states in the model multiplied by the total length of all the training sequences. Details are given in the excellent tutorial article onHMMs by Rabiner (1989). “he forward part of the forward-backward procedure can also be used to efficiently compute. -log Prob(sequencetmode1). the negativelogarithm of the probability of a sequence.given the model (as defined in equation (2)), without summing overall possible paths for the sequence (Rabiner, 1989). We call this the negative log likelihood (SLL)-score of the sequence. The average hTLL-acore of a training sequence is inversely related to the likelihood of the model, given by equation (3). and hence serves as a numerical measure of progress for each be iteration of the EM procedure. The KLL-score can also used to evaluate how well the model fits a novel “test.” sequence not present in the training set, as described in section (c) below. 1 ou5 mizing the probability of the path and then taking the negative logarithm. i t is convenient (and equivalent) to simply minimize the negative logarithmof the probability over all paths. This minimum we will call the distance from the sequence to the model, dist(s. model) = min {-log &ob(s, pathlmodel)) parh N+ = min. par& 1 2 [-log~(qiIqi-l)-log~(z,,,tg,)~ i-I This distance from a sequence to a model is analogous to thestandard“editdistance’’fromonesequence to another (nit,h gap penalties), see e.g. Waterman (1989), but isperhapsmore related tothedistance from a sequence to’a profile. The term -log9(rl~nlqi) represents a penalty for aligning the amino acid r,,, to the position represented by state pi in the model. The term -logF(qilqi-l) corresponds to a penalty for using the this is a transition from to pi in the model. If transition from a match state to a delete state, then this represents a gap-initiation penalty; if it is from a delete it representsagap-extension statetoadeletestate penalty; if it is from a match state to an insert state, it represents an insertion-initiation penalty; and if it is a transition from an insertion state to itself (a “self-loop”), then it. represents an insertion extension penalty. One of the main features of this distance measure isthat all t h e s e ’ penalties depend on the position in the model, whereas they would be fixed in most standard pairwise alignment methods. Often the most likely path has a significantly higher probability than all other paths, and in that case t.he distance defined here will be approximately equal to the SLL-score definedearlier. i Thecomputationtime for the Viterbialgorithm & proportional to the number of states in the model multiplied by the length of the sequence being aligned, i.e. the same as t.he time for the forward-backward algorithm. In addition, with a simple extension to the algorithm, the most. .probable path itself can be found using the usual backtrackingtechnique(Rabiner. 1989). This is the method we use to obtain our multiple alignments: each sequence is aligned to the model by the Viterbi algorithm, after which the mutua! alignment of the sequences among themselves is then determined.? (d) Using the H M M lo clwter sequenced and discover dfamilies When a relatively large numberof sequences are available, it is sometimes possible to obtain improved results (c) The Viferbictlgwithm and mzclliple alignmenf bydividingthesesequences into clustersofsimilar f r m an H M M sequences and training a differentHMM for each cluster/subfamily.The results of this are illustrated in more Theforward-backwardprocedureisrelated to the detail in Results section (a).Given a large set of unlabeled dynamic programming technique used to align one and nnaligned sequences,a simple extensionof the hidden sequence.to another, or more generally to align a sequence to a profile. A ,variantof the forward-backward procedure Markov model enables us t o use the EM training algorithm to automaticallypartitionthesequencesinto known as the Viterbi algorithm is similarto the standard By iterativelysplitting clusters of similarsequences. profile alignment algorithm (Waterman & Perlwitz, 1986; clusters, this method might be useful for building phyloBarton k Sternberg, 1990; Gribskov et al.. 1990). Instead genetic trees in a “top-down” manner. However, when of calculating t.he ELL-score for a sequence, which implithe clusters becometoo small there will be a n insufficient citly involvesall possiblepaths for that sequence through number of sequences in eachcluster to constructan the model, the Viterbi algorithm computes the negative logarithm of the probabilityof the single most likely path accuratemodel, so some”bottom-up”processingmay still be necessary. for the sequence. We can write this as In order to discover w clusters in the data, we make 2a - log max Prob(s, pathlmodel), (6) copies of the HMM, one for each cluster. We call these pork where Prob(s. psthlmodel) is given in equation (1). with u = z1...rLand path = qo ...qN+ l . Instead of first maxi- t We make no attempt to align portions of the sequence that use the insert statea of the model. 1506 zyxwvutsrqpo zyxwvu Hidden A1 czrkozy Models knownsubfamilies of the sequences. Esperimeotswith the clustering of globin sequences are described in Results section (a). zyxwvuts zyxwvutsrqponmlk nuu I (e) Modeling protein domains with an H l . t 1 nuu z BEGIN END There are many caseswhen one does not want to build a statistical model ofafamily ofwhole proteins like globins, but instead t o build a model of a structural motif or domain that oc*cursas a subsequence in many different kinds of proteins. such as the EF-hand motif (Kakayama et al., 1992) or the kinase catalyticdomain(Hanks & Quinn, 1991).Here we expect. our model'to only match a relatively small subsequence of any given protein, wit.h many other unmatched amino acids appearing before and after thissubsequence. One approacht o this problem is to alter the dynamic programming method used to align a sequence to a model so that. it tries all possible ways of aligningeachsubseyuenceofthesequence to a model (Waterman. 1989).lye use a simpler (but almost equivalent) method in which only the HMY model is altered. so thatthesamestandard procedures(forward-backward ant1 Viterhi) \vhich we use for models ofwhole proteins can he used \rithout modification for modets of domains. ('onsider a training set of many unaligned sequences consistingnot ofcompleteproteins.but of a specific domain.Our tirst step is to trainan HYM for these sequences esactly asdescribed earlier. As shown in Figure I. this HNJi \vi11 h a w initial and final "dummy" match states m, and mN- I (where S 1 = 5 in Fig. I ) that do not match any aminoacid. To alter theHJIM to represent a protein domain. we cntrrte 2 new insert states i, and i,. adding i, to the n~odelbefore the state nlo and i, at the end of'the model after mNTI(see Fig. 3). K P then acid a new clntnm_v BEOIS state before i, ant! a nrw duintny R S I ) state after i,. Eight new transitions are also atlclrtl to the IIIOCICI. The first 4 are from BEGTS t o i,. fronl m y to i,. and the self-loops from iBto itself and froln i, to itself: Thew all haw the same prohabi1it.y 1'. for sonw 1' lwtwwn 0 and I . The second 4 transitions are fron~i, t o m,. from I3E:CTS to m,. from i, to ESD. and from 111,~- t o KXD. Thme a11 haw probability I p . The new states atltletl hrforr and after the modet. along with these transitions. form 2 new tnddules. 1 for matching theestra atnino acids that occur in the sequenc*..I)efi)re the hnain. and the other for matching the amino acids aftrr the clotnain. The choice of the Imranleter p does affrct the way that the overall rnwlel aligns with a given q u e n c e . To see ho\v. it is conwnirnt t o think of the negative logarithm of the prol)al)ility of a transition as a penalty for using that transition. as tlescrilwtl i n section (c). above. In the moditied n1cwlel. all sequences must suffer a penalty of -log ( I - p ) to enter and again to esit the domain of part the model. no matter which path they take. Hence this penalty is a tisrtl c w t . which ran fw ignoredwhen zyxwvutsrqpon zyxwvuts nuu Figure 2. HMJI architecturefordiscoveringsub- zyxwvutsrqp families. components of the(composite) HYJI. Presently. the number IC of clusters and the initial lengthsof the ~notlels €or these dustem are determinedempirically. \Ye then add 8 new begin state with w outgoing transitions. one t o each of the begin states of the component HWls (sre Fig. 2). This new begin state is analogous to the other begin states in that it generates no amino acid. IVe then train this composite model with the El1 algorithm as c1eucril)ecl in section (b). above. The EM re-estimation of a co111ponent model is the same as the reetimationof a sinplr model, except that the weight that a sequence has in the re-estimation of a component is proportional to thep r o b ability of the sequencegiven that com1x)nent nwclel. \ Thus, sequence that have betterS L I ~ c o r e for s a partih l a r HMN component hare greater influence in re-(?;timating the parameters of that component.and this ~ a u . w s t h e parameters of that component to change in such a way that thecomponent further *'specializes" in n~twleling those sequences. The "surgery" I)rowclurr clrscribetl belor in section (8) is used to adapt the length of that component to further specializeit. In this manner. the individual romponents evolve during training t o wprrsrnt clusters in the trainingsequences. This wayof using EN is called mixture modeling in the statistics literature(1)ucla gL Hart. 19f3: Ereritt Pr Hand. 1981). and is known as "(soft) competitive learning" in t h e neural nrtuork literature (Sowlan. 1990). When the model is trained.the probability of a sequence gi\-en any of the submodels van be ralruiatrtl. Le. the probability thatthesequenc! belongs tothe corresponding cluster/subclass. The negative logarithm o f this probability corresponds to the SI,1,-score calculatrcl for a simple HJI3I. As with the standard H W I weuse. this yields a quantitative measure of how \vel1 the n~odel fits the data. The dusters found can also lw compare<l to zyxwvu + -, Figure 3. H U M architecture for modeling donlains. - zy zyxwvutsr zyxwvutsrq zyxwvutsrq zyxwvutsrq zyxwvutsrqpo zyxwvutsrq zyxwvutsrqpon zyxwvut comparing the distances or S L 1 ~ c o r . eof~ 2 sequences with respect to the model. In addition to this penalty. all sequences willsuffer a penalty of /i ( - l o g p + log 20). where K 2 0 is the number of amino acids that are not matched to the original domain model.but are instead matched in the states i, and i,. The -log 20 term arises because we set t.he probabilities of each amino aridto 1/20 in the insertion states i, and i, (see Krogh el nl.. 199.3~). Thus p wilt determine the "pressure" on the sequence to align something to the domain model. i.e. if 1) is low it is advantageous to squeeze manyaminoacidsintothe domain model. using the insert states in this part of the model. If it is high. it is possible that most sequences in the wouldprefer to pass throughthedeletestates domain model. aligningex-erythinginstead to t h e new modules before andafterit. It is straightforwardto estimate p the same way a s all the other paranwters. the only additional problem is that the rrilnw value nmst he used in all the transitions that 11.w this v a l w . -tying" model theseparameters to eachother.Othrnviscthe might become biased towards aligning the tlomain either near the beginning of the sequence or near the end of the sequence. We hare not attempted t o estimate p. Rather. we have useda fixed p = 1 with goor1 results. (This should be thought of as a limit of p approad~ing1. otherwise -log (1 p ) is infinite.) Usingthisconstruction. itmay also he 1)ossible to discover interesting domains by training on whole protein sequences, and letting EM determine which part of the proteins to model. Furthermore. if more than one occurrence of t.he same domain is expected in sonle sequences. then this model can be further modified to find all occurrences. This is accomplished by simply addinga transition from the EXD state back to the.BE(;IS state. - (f) Searching a database wifh an H.11 SI Once an HMM is built for a family of proteins. it can be used to search a database such as PIR or S\VISS-PROT for ot.her proteins in this family. Similarly.if an H l l l i is built for a protein domain or motif. then it can be used to search for occurrences of this domain or motifin the database. much like a PROSITE espresion (Rairoch. 1992). a commonly used method for searching for patterns found in protein sequences. Like a profile ( f a t e m a n (I: Perlwitz. 1986 Barton & Sternbeg. 1990: Gribskov el a/.. 1990: Bowie et ai.. 1991: Liithy et al.. 1991). an HJIJI has an advantage over a PROSITE expression for database searches. It takes into account a large amount of statistical information in matching a sequenw. and weighs this informationappropriately.ratherthan relying on relatively rigid matching rules. As described in section (b). above. the fonvard part of theforward-backward'dynamirprogrammingmethod calculates a BLtscore for anytestsequencethat measures how well it fits t h e model. This SLL-scoreis the negativelogarithm of t h e probability of the sequence given the model. It turns out that this raw SLL-score is too dependent on the length of the test sequence to be used directly to decide if the sequence is in the family modeled by the HMIM or not. However. we canovercome this problem by normalizing this SLL-score appropriately. Whenever we build an HMM for a family of proteins or for a protein domain. we run all theproteins in a standard database (forinstance. SWISS-PROT) through thisHMM and compute the SLL-score for eachsequence. A scatter plot of sequence length r w w s SL1.-score for our kinase catalytic domain model is given in Figure 9. Most proteinstend tn lie on a fairlystraightline (towards the top of the plot) indicating that the is proportional to their SLL-score for theseproteins lengths. These proteins are the ones that do not contain the kinase catalytic domain and thus look like "random proteins" tothe kinasecatalyticdomain model. In contrast. the proteins that do contain the kinase catalvtic domaintend to have SLL-scores that are muchlower than espected for. proteins of their length. and hence appear below the linear band of non-kinase proteins. We can quantify thedifference between SLL-scores for prot.einscontaining the kinasecatalyticdomainand ST,I,-scores for proteins not containing the domain by a simple statisticalmethod, as follows. Csing a local windowing technique.t we first calculate a smooth average curve for the roughly linear band of the SLL-score rer.ms lengthplot. Thestandarddeviation around this average curve is also calculated. Using this. \vc. calculate the difference between the SLL-score of a sequence and the average SLGscoreof typical sequences of that ~ ~ length. n emeasured in standarddeviations. This number. is called t.he 2-score for the sequence. We then rhoow a 2-score rut-off. either [I priori or by looking at thehistogram of 2-scores for sequences in t.he database (see Fig. IO). and use it to decide if a given sequence fits the model or not. We havefoundthat a 2-score of approsimately .iappears a good choice in most cases we hare esaminrtl. but we suggest carefully checking the histogram by eye before deriding oncut-off. a For esample. for our HMM of the kinase catalytic domain. sequenceswith 2-scores below 5 are classified as not 'c-ontaining the kinasecatalyticdomain,and sequences 5 are classified as containing $e with2-scoresabove catalytic domain. If the 2-score of a sequence indica& that it contains the catalytic domain. we can align the sequence to the catalytic domain HMlI to find out where this domain occurs in the sequence. The time it takes to do a database search is proportional to the number of residues in the database times the length of the model. For our globinmodel(length 145) we cansearch the 8.375.000 residues)in WISM-PROTdatabase(about approsimately 2 CPZ: hourson a Sun Sparcstat.ion 1. I'sing the shorter EF-handmodel (length 29) it takes only 18 C'PV .seconds (1 1 user min) on a Sun Sparcstation2. X parallel implementation of the search procedure (not r e t implemented) r i l l speed up these searches substantially, as it has the EM training procedure. While the statistical techniqueswe have used to determine 2-scoresare still quite crude.we hare found that the HYJIs are sufficiently good models that these techniques work \vel1 enough in practice. However. i t may be that more sophisticated techniquesare needed in certain cases. ' zyx t The arerage curve is calculated a s follows. For each length i starting at i = 1. the length I, is computed such that there are at least 500 proteins of lengths i to li and less than 5013 proteins of lengths i to I, 1. The length interval i to li is called a window. The average curve is piecewise linear through the pointscorresponding to the average length and average XU-score for each window. The first and last parts of the curve are calculated by linear regression in the first and last window. respeetirely. The standard deviat.ion of t.he points from the smooth curve is also calculated for earh window. The estimate of t.he average run.e can he impmwd by eliminating outliers. i.e. XLT.-scores that lie many standard deviations from the average. W e iterate the p r o w s of removing outliers and re-estimating the average curve until no more outliers remain. - zyxwvutsrqp zyxwvutsrqpo zyxwvuts zyxwvuts zyxwvu zyxwvu zyxwvut zyxwvutsr zyxw (g) Initial d e l , i d minima, and choice of model length As mentioned in section (b), above,when estimating the model from the training sequences, the EM algorithm does not guarantee convergence to the best model. It is basically a steepestdescent-type algorithm that climbs the nearest peak (local maximum) of the likelihood function (or the posterior probability in MAP estimation). Since finding the globally optimal model seems to be a difficult optimization problem general in (Abe & Warmuth, 1990), we haveexperimentedwithvarious heuristicmethods to improve theperformanceofthe method. Probably the best method is to give the model a hint if something is already known aboutthe sequences, which is often the case.A good starting point makes it much more likely that the nearest peak is at least close to optimal. This is donebysettingtheprobabilitiesintheinitial rnodei to values reflecting that knowledge. If, for instance, an a l i m e n t of some of the sequences is available, i t is straightforward to translate that into a model by simply cnlcuIating the relative frequencyof the amino acids and the transitionfrequenciesineachposition, as in the profile method (Gribskov et al.,1990). It is of course even more interestingif the model can be found from a M a ‘ m a ,i.e. using no knowledge about the sequences. For that we haveusedaninitialmodel where allequivalentprobabilities are the same,i.e. d ( q + l I m k )is independent of the position k in the model, and similarly for all other transition probabilities, and 9(zlmk) is also independentof k. To avoid the smaller local maxima. noise is added to the model during the iteration before each re-estimation. Initially quite a lot of noise is added, but over 10 iterations the noise is \decreased linearly to zero. Since noise is added directlyto ‘$he model, it is not like the usualimplementationof iimulated annealing, but the principle is the same. The “annealing schedule” is p&ntly rather arbitrary, but it does seem to givereasonableresultsf if it is applied several times,and thebest of the models found is used as the final model. It is important that the best model be selected, since suboptimalmodelsdoproduceinferioralignmentsin general. However, when studying alignments from suboptimal globin models, we noted that they tend t o align some regions well. occasionally getting better alignments in those regions than the best overall model found, while in other regions they are completely incorrect. This leaves open the intriguing possibility of combining the best solutions found for diRerent regions into a new overall best model. We have not yet explored this possibility. The lengthof the model is also a crucial parameterthat needs to be chosen a priori. However, we have developed a simpie heuristic that selects a good model length, and even helps in the problem of local maxima. The heuristic is this: after learning, if more than a fraction$ ydel of the paths of the sequenceschoose dk, thedelete state at the model. position E , that positionisremovedfrom Similarly, if more than a fraction yinr make insertions at. position k (in state it), a number of new positions equalto the average number of insertionwmade a t that position are inserted into the model after position k. After these changesin the model, i t is retrained, and this cycleis repeated until no more changes are needed. We call this “model surgery”. (h) Over-@ing and H A P mtimatia A model with too many free parameten cannot estimated well from a relatively smalldata set of training sequences. If we try to estimate such a model, we run into the problem of overfitting, inwhich the model fib the trainingsequencesvery weI1, but gives a poor fit to related (test) sequences that were not included in the training set.We say that themodel does not “generali~” well to test sequences. This phenomenon has been.well documented in statistics and machine learning (see e.g. Geman el al., 1992; Berger,.1985). One way to deal with this problem is to control the effectivenumberof free parameters in the model by using prior information. This canbeaccomplishedwith MAP estimation. Parametem that we assume (viaour prior distribution onmodels) ~ 8 n be wellestimated a priori in effect become less adaptive, because.it takesa lot of data to override OUF prior beliefs about them,whereasthose about whichwehaveonly weak prior knowledge are estimated in almost the same manner as inmaximumlikelihood estimation. In this way. the model can have a very large number of para. meters, but a much smalier number of “effectively free” parameters. To make MAP estimation practical, we use Dirichlet distributionsas priors. Thedetails of the method are described elsewhere (Kroghet al., 1993a; Brown d d., 1993). 3. Results (a) Globin experiments The modeling was first tested on the globins, a large familyof heme-containing proteins involved in the storage and transport of oxygen that have different oligomeric states and overall architecture (for a review see Dickerson & Geis (1983)). Hemoglobins are tetramers composed of two OL chains and two other subunits (usuallyB, y , d or 8). Myoglobin is a single chain, some insect. globins are present as dimers and someintracellularinvertebrate globins occurinlarge complexes of many subunits. Globin sequences were extracted from the SWISS-PROT database (release 19) by searching forthekeyword “globin”. Eliminating the false positives. resulted in 625 genuine globin sequences of averagelength 145 amino acids. We left three non-globins in the sample for illustrational purposes giving a total of 628 sequences. Thesample of globins in the database is not the random sample a statistician would prefer, but is perhaps one of the bestandlargest collections of protein sequences from a homologous family. Sesrching for the words “alpha“, “beta”, L‘gamma’’, “delta”, “theta”, and “myoglobin” in the data file yielded 224 alpha, 199 beta, 16 gamma, 8 delta and 5 theta chains and 79 myoglobins, which adds up to 531 sequences. These should naturally be considered minimum numbers. but they give a good picture of howskewed the sample is. To t e s t our method, we trained an H M M using the method described in Methods sections (b) and zyxwvutsrqp t An alternate method that also appears to give good results has been developed by Baldi el al. (Baldi el al.. 1993: Baldi & Chauvin, 1993). This method uses stochastic gradient descentin place of the EM method, which may help in avoiding local minima. 1Currently we choose yde, and yinr each to be 1/2. zyxwvutsr zyxwvutsr zyxwvu zyxwv zyx zyxwvuts zyxwvutsrqp H zdden ill arkoa 111 odels 1DUY Helix AAAAAAAAAAAAAAAA BBBBBBBBBBBBBBBBCCCCCCCCCCC DDDDDDDEE HBA-HUNAN ---------V L S P A D K T N V K A A W G K V G A - - H A G E Y G A E A ~ ~ F ~ F P T - - -BGSA -HBB-HUNAH --------VHLTPEEKSAVTALWGKV----NVDEVGGEALGRLLWYPWTqRFF~FGD~TPDAVHGNP MYG-PHYCA ---------V L S E G E W q L V L H V W A K V E A - - D V A G H G q D I L I R L F K S H P E ~ ~ F D ~ H ~ ~ A E H K A S E GLB3-CHITP ----------LSADQISTVqASFDKVXG------ DPVGILYAVFKADPS1HP;K.FTQFAG-KDLESIKGTA GLB5,PETMA PIVD'PGSVAPLSAAEK?XIRSAUAPVYS--nETSGVDILVK~STPAAqE~PK~GL~ADQLKKSA LGB2-LUPLU --------GALTESqAALVKSSWEEFNA--NIPKHTBRFFILVLEIAPAAKDLFS-nK-GTSEVPqNNP GLBI-GLYDI --------- GLSAAQRqVIAATWKDIAGADNGAGVGKDCLIKFLSAHP~HAAYFG-FSG----AS---DP EEEEEEEEEEEEEEEEEEE Helix FFFFFFFFFFFF FFGGGGGGGGGGGGGGGGGGG HBA-BUNAB QVI(GBGKKVADALTNAVAHV---D--D#PNALSALSDLHAHKL--RVDP~YFKUSBCUVTLAAELP~ HBB-BUMAH K~AH~KKVLGAFSDGLAHL---D--NLKGTFATLSELHCDKL--HVDPEHFRLLGBVLVCV~AEEFGKE HYG-PHYCA D L K K B G V m L T A L G A I L K K - - - - K - G H H E A ~ P ~ q S H A T K H - - K I P ~ ~ ~ I S ~ I I ~ V L H S ~ ~ P G D GLB3,CHITP PFETBANRIVGFFSKIIGEL--P---NIEADVNrrVASHKPRG---VTHDq~N~~AGFVS~HT--D GLB5,PETMA DVRWBAERIINAVNDAVASM--DDTMHSIIKLRDLSGKHAKSF--q~P~~~AVI~TV~G---LGB2,LUPLU ELqAaAGmr~LVYEAAIqLQVTGVWTDATLKNLGSVHVSKG---VAD~P~EAILKTIK~VVG~ GLBI-GLYDI GVAALGAKVLAqIGVAVSHL--GDEGKllVAqHKAVGVRHKGYGNKHIKAQYFEPLGASLLSAHEHRIGGK ' ' Helix HHHHHHHHHHHHHHHHHHHHHHHHHH HEI-EUKAI FTPAV3ASLDKFLASVSTVLTSKYR-----BBB-HUZIAB FTPPVqAAYQKWAGVANALAHKYH-----HYG-PHYCA FGADAQGAHHKALELFRKDIAAKYKELGYQG GLB3,CHITP FA-GAEAAUGATLDTFFGHIFSKM------GLB5,PETMA DAGFEKLnSHICILLRSAY------LGBZLUPLU WSEELlDSAVTIAYDELAIVIKKEMNDAA--CUI-GLYDI MUAAAKDAUAAAYADISGALISGLQS----\ Figure 4. Seven representative globin sequences of known structure and their alignment taken fmm Bashford et a.S, ----- zyxwvutsrq (1987). The letters A to H in Helix denote the 8 different a-helices. Some regions, especially CD,D and FG, are not well defined. The sequences and their SWISS-PROT identifiers are Human a (HBA-HUMAK), human fl (HBB-HUMAN), sperm whalemyoglobin (MYGSHYCA), larval Chir0nomou.s thvmmi globin(GLB3-CHITP). sea lamprey globin (GLWETMA), Lupinus Z u & m leghemoglobin (LGBLLUPLV): and bloodworm globin (GLBLGLYDI). (In SWISS-PROT 19 a S is used instead of an ''-" in the identifiers.) (gf. We usedahomogeneousinitial model that contained noknowledge about theglobin family. Its probability parameterswere derived from the prior, and were the same for all equivalent transitions (i.e. 9 different transition probabilities). All amino acid probabilities (theB distributions) were set equal to the distribution of the amino acids given by Krogh et d . (1993a). In the insert states we used a probability of 1/20 for all amino acids. The only model parameters set by handare the initial transition probabilities and correspondingregularization parameters (see Krogh et d.,1993~).From our experience, the method does notseem to be very sensitive to the choice of theseparameters,but i t would requireconsiderable furtherexperimentationto verify this quantitatively. For our training set, we picked 400 sequences at random from the 628 sequences. We withheld the remaining 228 sequences in order t o test the model on data notused in the training process Themodel was trained using noise and model surgery (ysrr = y h = 0 5 ) , as described in Methods section (g). This procedure was repeated about 20 times with model lengths chosen r a n d o d y between 145 and 170. The average run-time was around 60 CPU minutes on a Sun Sparcstation I. For each run we computed a XLL-score for the model, which was the average of for thetraining sequences, as theELL-scores defined in Methods section (b).The final NLL-scores varied considerably for these runs but the best was 210.7. We then took this model, produced ten new models by adding noise, and optimized these. These models all generated approximately the same NLL-score and we picked the model with the best NLL-score, 2103, having a length of 147. W e validated this modelt in two ways: from the alignments it produced, andbyitsabilitytodiscriminate betweenglobins and non-globins.The results are described below. (i) Multiple sequence a l ~ m e n t s A multiple alignment of many globin sequences has been produced by Bashford et d. (1987) by including intothealignment procedure tertiarystructure information of seven globins (Fig.4). This t We stress t h a t the final model was chosen according measure, namely the KLL-score on the training set. and not retroactively on the basis of how well it did in multiple alignment or database search to an objective tasks. zyxwvutsrqpo zyxwvutsrqpon zyxwvutsrq zyxwvut zyx zy 1610 Markov was achieved by aligning these seven sequences and then aligning the rest of the 226 studied to the closest of these seven. Incontrast. generating multiplealignments with HMMs requires no prior kno.rvledgeof underlying structure. Using the globin HM3fCI, we produced a multiple alignment of all the 625 globin sequences by the Viterbialgorithm as described in Methods section (c).Figure 5 shows this alignment for the seven sequences from Bashford et a l . (198i). HBA-HUHAN HBB-BUXAN MYG-PHYCA GLBI-CHITP GLBS-PETMA LGB2,LUPLU GLBI-GLYDI \ '\. AAAAAAAAAAAAAAAA F- HBA-BUHAN HB3,HLJHAN HYG-PHYCA GLB3,CHITP GLB5-PETHA LGB2,LUPLU GLB1,CLYDI BBBBBBBBBBBBBBBBCCCCCCCCCCC DDDDDDDEE + V.........LSPADKTNVKAAWGKVGA..HAGEYGAEALERHFLSFPTTKTYFPHF-DLSHGSAQ---Vh........LTPEEKSAVTALWGKV--..MmEVGGEALGRLLWYPWTQRFFESFCDLSTPDAVMGNP V.........LSEGEWQLVLAVWAKVEA..DVAGHGqDILIRLFKSBP€TLEKFDRFKHLKTEA€~ASE -..........LSADqISmQASFDKV--..KGDPVG--ILYAVFKADPSIMAK~~-AGKDLESIKGTA PivdtgsvapLSAAEKTKIRSAUAPWS..TYETSCVDItVKFfiSTPAAQ€~~PKFKGLTTADqLKKSA Ga ........LTESQAALVKSSWEEFNA..NIPKHTHRFFILVLEIAPAdKDLF-SFLKCTSEVPq-NNP G. . . .LSAAq~VIAATWKDIAGadNCAGVGKDCLIKFLSAHP~MAAVF-GF----SGASD--'P . . . .. FFFFFFFFF ********* +****************** Helix (ii) Database search.: discriminating globins f.Om non-globins The globin HMM model we found was also tested on ail t.he 25.044 proteins in the SWISS-PROT database release 22-0 of length less than So00 amino acids (which is all but 2). A XLL-score and a 2-score were computed for each of these sequences as described in Methods section ( f ) . These are plotted in Figures 6 and 7 as a scatter plot and a histogram, respectively. For the histogram (but not the scatter plot), the data were filtered as follows: All sequences witha2-score >3.5 andeither more than a total of 23, or more than IS*/,unknown residues were removed (a total of 23).Currently, we treat an unknown a.mino acid, X,as being the most probable amino acidat theposition it is matched to, ***************+ ++++++++****************++* EEEEEEEEEEEEEEEEEEE Helix HBA-IIUHAI HBB-EUHAN WG-PEYCA GLB3-CHITP GLB5-PETHA LGBP-LUPLU GLBI-GLYDI between secondary structure elements, The 1 s t two insertions appear in the F/G region. zyxwvuts zyxwvutsrq The alignment found in ,this experiment agrees estremely well with the structurally derived alignment of Bashford et al. Our alignment differs in the region between the C and E helices. However, this is a highly variable areasinceonly some globins possess a D helix. The difference in the F/C-helices is mare pronounced. with the remaining discrepancies possibly representing an alternative alignment. Four of the insertions the model chose are in variable regions between or at the end of helices.i.e. Helix Modds Hidden FFFFFGGG GGCGGGGGGGGGGGGG ***************+ -VKGHGKKVADALTNAVAHVDD.....HPIALSALSDLHA...HKLRVDPV.NFKLLSHCLLVTLAAHLP KVKAHGKKVLGAFSDGLAHLN... ..LKGTFATLSELHC ...DKLHVDPE.NFRLLGN~VCVLABHFG DLKKHGtlTVLTALGAILKKKGH.....HEAELKPLAQSHA...TK-HKIPILYLEFISEAIIHVLHSRRP PFETHAPRIVGFFSKIIGELPN . .IEADVNTFVASHK.. .PR-GVTHD. QLHIFRAGFVSYMKAH-DVRWHAERIIIAVNDAVAS~Dtek..HSXKLRDLSGKHA ...KSFqVDPq.YFKVLAAVIADTVAA--ELPAHACKVFKLVYEAAIQLQVtgvvvTDATLKHLGSVKV...SK-GVADA.HFPWKEAILKTIKEVVG GVAALGAKVLAQIGVAVSHLGDegk..MVAqPIXAVGVRHKgygNK-HIKAq.YFEPLGASLLSAMEHRIG ... HHHHHHHHHHHHBHHHBAHHHHHHHH ++**+****************** AEFTPAVHASLDKFLASVSTTLTSKY......R KEFTPPVQAAYQKVVAGVANALAHKY ...... H GDFGADAQGIMlKALELFRKDIAAKYkelgyqG TDF-AGAEAAWGATLDTFFGMIFSKH...... GD------ AGFEKLMSMICILLRSAY. ..... AKWSEELRSAWTIAYDELAIVIKKEMnda ...A GKHNAAAKDAWAAAYADISGALISGLq .....S - - Figure 5. The alignment Of the same 5 globins as in Fig. 4. as obtained from our model trained on 400 randomly chosen globin sequences. The capit.al letters represent amino acids aligned to themain line of the Inoclel. -. t o deletions in the model. and lower-case letters to amino acids treated as insertions by the model. The . is used as a fill character to accommodate insertions. S o attempt has been made to align the insertion regions. Tn the line above the alignments * indicates complete agreement of a column ivith the structural alignment (Fig. 4) and + denotes a minor deviation (the only accepted difference is a reasonable displacwnent of a gap). The regions between the helires are not checked in this way. The training set contained 5 of the 5 globins. not HBX-HUMAS and GLRLPETMA. zyxwvutsrq zyxwvutsr zyxwvut zyxw zyxwvut zyxwvutsrqpo zyxwvutsrq l 100 1 . , , Length of Sequence I 200 , ( . . , 300 Figure 6. Plot of SLL-SCOW. wr.wx sequentv length for globins an? non-globins. ..MI sequences of length less than 300 from the S\YISR-PROT P i database are sho\va. including partial sequences and3 false globins from the globin file. and sequences containing database from the man?- Ss. \ E3 non-globins training set Itest set 3000 10 e 6 4 2 0 Z-score Figure 7. Histogram showing t h e number of sequences with a certain 2-score. The trainingset of 397 globins. the test set of 231 globins. and therest of thesequences from SWTS+”-PROTP2 after “filtering” are shown. The insert shows expansion of the region around a Z-score of 6. so sequences with many lis spuriousi~match the model very well. release of Since we searched a newer SWISS-PROT (release22) than the one from which the globin training set was extmcted (release 19), eight new globins were found and incorporated into the test set. Five globin fragments of length 19 to 45 were removed from t.he data. Three non-globin sequences in the globin file that were identified as outliers in Figure 6 were removed. One of thesenon-globins wasleft as part of the trainingset toillustratetherobustness of the method. The model distinguishes extremely well between globins and non-globins. Choosing a 2-score cutoff of 5 we would miss 2 out of 628 globinst and get essential]?nofalsepositiveglobins.Thereis one ”non-globin”, a bacterial hemoglobin-like protein (SWISS-PROT id HMP-ECOLI)l that may or may not be-countedasa false positive. Only one sequence, the heme containing catalase of zyxwvu zyx Penicillium vilale (CATA-PENVI, Z-score C7),has a Z-score between 4.2 and 5-1,so any cutoff in this range will essentially give the same separation. The two sequences falling between a. Z-score o f 1. and 4 t 628 in the original data set. plus 8 new. minus 3 spurious. minus b fragments. 39i were left from t h e training and the remaining 231 made up the test set. zyxwvutsrqpo zyxwv zyxwvutsrq zyxwvuts zyxwvutsrq zyxwvut 1512 Hidden Markov Models (GLB-PARCA and GLB-TETPY) are protozoan, whereas the other globins are metazoan. The primary sequences of these globins are similar and have little similarity with othereukaryotic globins. Note also that both of these sequences are in the test set. (iii) D i m e r i n g subfamilies of globins We also performed an experiment to automatically discoversubfamilies of globins using the method described in Methods section (d). An HMM with ten component HMMs was used. The initial lengths of the components were chosen randomly between 120 and 170, but were adjusted by model surgery during training. We trained this HMM on all 628 globins and then calculated the NLL-score for each sequenceforeach of theten component HMMs. A sequence was classified as belonging to thecluster represented by the component HMM that, gave the lowest KLL-score, i.e. the one giving the highest probability to t h a t sequence.? Three of these clusters were empty and the remaining seven nonempty ones represented chains from known globin subfamilies: C b s X. 233 sequences: principally all CY,a few [ (an a-type chain of mammalian embryanic hemoglobin), %/x' (the counterpartof the a chain in major early embryonichemoglobin P), and 8-1 chains (early erythrocyte a-like). Class 2. 232 sequences,almost all a few 6 \(p-like), E (&typefoundinearly embryos), y (comprise fetal hemoglobin F in combination with 2 Or chains), p (major early embryonic /.?-type chain) and 6 chains (embryonic &type chain). C b s 3. 'il myoglobins. C h s 4. 58 sequences. The 13 highest scoring in this cluster are leghemoglobins. This class contains a variety of sequences including the three non-globins in original data set. Class 5. 19 sequences. Midge globins. Class 6. Eight sequences. Globins from agnatha (jawless fish). Class 7. Seven sequences. Varied. We havenotrepeatedthisexperiment using different randomization t o ascertain if better results can be obtained. However, we are encouraged by the results of this first experiment since it is able to classify correctly t h e major globin subfamilies (alpha, beta and myoglobin). indicates what fraction of the 400 training sequences made that transition or used that padicular amino acid. A broken line indicates that less than 5% of the sequences used that transition. (The continued delete is mostly due to fragments that have to make many deletions.) The histogram in a match state shows the distribution of amino acids that were matched to that state. The number ina n insert shows theaverage length of a n insertion beginning at that position. Forthe amino acids the orderingproposed by Taylor (1986) is used. S i r t i n g from the top, the amino acids are medium-sized and non-polar, small and medium polar (around G and P), medium sized and polar (around K), large medium-polar (around F and Y):and finally below they are medium-large and non-polar. There does seem to be some tendency for t.he distributions to peak around neighboring amino acids when using t.his ordering,.as one would expect. When one looks at the whole model, regions that are highly conserved are also readily distinguished from the more variable regions, both as a function of the probability that. a position is the entropy of thedistribution of skipped,and amino acids at that position. zyxwvu zyxwvu . (iv) The $rial globin model Examination of the model itself yields information on thestructure of globins.Figure 8 shows t h e normalized frequency counts (the numbers used to re-estimate the parametersof the model) from some parts of the final model. The thickness of a line ~~ 7 We can also calculate the posterior probability of each cluster by looking at the transition probabilities out of the global start state, and thereby obtaining a posterior distribution over the 10 clusters for each sequence. However, these posteriors are very sharply peaked, so this adds little to t h e analysis. (b) Kinase experiments Protein kinases are defined as enzymes that transfer a phosphate groupfrom a phosphate donor onto an acceptor amino acid in a substrate protein (Hunter, 1991; Hanks et d.,1988). Based upon the accept.or aminoacid specificity, theyhave been classified intoserine/threonine,tyrosine,histidine, cysteine, aspartyl and glutamyl kinases. Only enzymes in the first two categories have been well characterized and recent deyelopments indicate that some can phosphorylate both a.lcoho1 (serine/ threonine)and phenol (tyrosine)groups, the socalled dual-specificity protein kinases (Lindbmg et a.1.: 1992). It. is the region comprising the catalytic domain of thesehydroxyaminoacidphosphorylating enzymes t.hat we model by an H M M and which we subsequently refer to as protein kinases or simply kinases. Despite the differences in size, substmte specificity. mechanism of activation, subunit composition a.nd subcellular localization:all these kinases share a homologous catalytic core c,ontaining 12 conserved subdomains or regions (Hanks 8; Quinn, 1991; Hanks et al., 1988). Because the kinase catalyticdomain is onlya subsequence embedded in a Iarger protein,the kinase experiments differed from the globin experiments. The HMM used in the globin experiments modeled t.he entireprotein ratherthansimply a segment of a protein as is the case for the kinase family. Modeling domains requires several modifications t o ourstandard HMM training which are described in Methods section (e). The training set. for these experiments is a group of 193 sequences from the March 1992 release of the protein kinase catalyticdomaindatabase maint.ained by Hanks & Quinn (1991). Thisset is Hidden Markov Models zyxwvu zy 1513 zyxwvutsrqp zyxwv zyxwvu zyxwvuts Figure 8. Parts of the final globin model. The position numbers are shown in the delete states. ' composed of serine/threonine,tyrosineanddualspecificity kinases principally from vertebrates and higher eukaryotes but also includes some from lower eukaryotes and viruses. We trainedten HMMs on all 193 (unaligned) sequences in this-data. set using the prior distributions described by Krogh et ul. (1993~).N o parameters. of the modeling process were set manually and the initiaI model lengths rangedfrom 242 to 282 positions (this encompasses the average length of the sequences in our kinase catalytic domain training set). At the endof the ten training runs, the best kinase model had a NLL-score (theaverage --logP(sequenceImodel) overthetrainingset) of 588.39 and a length of 254. Modules were added at the beginning and end of this model as described in Methods section (e). We tested this modelin the samemanner as described earlier for the globin model. Our maintests were discrimhation tests,in which we utilized the model. to search the SWISS-PR.OT version 22 database (25,044 sequences) for proteins containing the kinase catalyt.ic domain. As described in Methods section ( f ) , a NLL-score was computed for each of the sequences in the database and this information was used to compute a sequence's deviation from the average curve as measured by a Z-score. The data were then filtered to remove all sequences with any unknown residues (353) and all sequences having length less than 200 (4230): since complete protein kinase catalytic domains range from 250 to 300 residues (Hanks et al., 1988). This filtering removed a total of 4386 sequences. A scatter plotof NLGscore v e m w length for the SWISS-PROT sequences is given in Figure 9. A cutoff of 6-0 was chosen because there are no sequences with Z-scores between 4935 and 6773. See Figure 10 for a histogram of the resulting Z-scores. Any sequence having a Z-score >6.0 was therefore classified as containing the kinase cata.lytic domain while thosewith Z-score t6.0 were classified as not possessing the domain. With this cutoff, 296 sequences were classified a s containing the kinase cat.alytir domain. The remaining 20,357 srquerx-w \\-ere rejecttd. 1514 zyxwvut zyxwvutsrqpon zyxwvu Hidden Yarkoa Models t zyxwvutsr zyxwvuts , Certain Kinases . Rest of SWISS-PROT zyxwvutsr ~ ~ ' " ' 500 ' ' " " " 1000 ' ' " ' t1500 2000 Length of Sequence Figure 9. scatter plot of SLL-score \-ersus length for sequences in SiVISS-I'ROT using the Kinase HMX ', The genera.1 issue of estimating the number of false negatives and false positives when distinguishing sequences belonging to a given family ICertainKinases a * . All other sequences in SWISS-PROT . from non-members is a complex one. In the case of 'the glohins. it is "relatively" straightforward since it 4s possible to identify all the globins in t.he database l ) prfonning ~ a keyword or title stringsearch. The situation for the kinase domain or the EF-hand motif (see section (c) below) is less obvious and thus more problematic. For instance, while a given protein may possess the' sequencecharacteristics for this motif or domain. functionally, the region may not bindcalcium or possess kinaseactivity. We hare attempted to address this complicated matter as best we can as described below. However, we stress that ~ r - edo not feel able to give a definitive answer as to the number of true false negatives and true false positives in our kinase or EF-hand database discrimination test.s. A listof potentialprotein kinases was created from the unionof sequencesdesignatedas being kinasesfrom four independent sources: our HMM, PR.OSTTE (a dictionary of sitesand patterns in proteins(Hairoch. 1992)), PROFTLESEARCH (a technique used to search for relationships between a protein sequence and multiply alignedsequences (Gribskor et ai., 1990)) and a. keyword search. Two regions of the catalyticdomain of eukaryotic protein kinases have beenusedt,obuild PROSITE signature patterns. The first pattern corresponds to an area believed to be involved in ATP binding (PROSTTE: entry PROTEIXKIXASEATP, sequence motif [ l,I\-]G.C.tF~~lJ[Sc].\'). There are two signature patterns for the second region important for ratalytic activity: one specific for serine/ zyxwvu ClO 0 10 20 30 40 ! Z-score Figure 10. Histogram showing the number of sequences with H certain %-score relative to the kinase model. Hidden Markov Models zyxwv zy 1515 zyxwvutsrq zyxwvu zyxwvuts OK BOG). bovine cGMP-dependent protein kinase kinases (PROTEIS-KTSASE-ST, (OKB02C'). bovine protein kinase C (KIBOC), ILIV~~FYC].[HY?.D[LIV,MFT]~.PS[LIV.MFYC]I) threonine human mos kinase-related transforming protein and the other for tyrosinekinases (PROTEIS(TVHVFG): human r e f a kinase-related transKIPU'ASE-TTR, [LIV~F~C}.[HY].D[L~~~lF~J [RSTA].ZX[LIVXFYCJB). Since PROSJTE formingprotein (TVHK"), mouse pim-I kinaserelated transforming protein(WMSPI),and human expressions do not albw for flexiblegapping or kinase-related transforming protein insertions, a profile of k i n a was constructed from feslfps (TVHVFF). The keyword search consisted of an alignment of seven kinasesand employed for database discrimination tests (11. Gribskov, sea.rching t.he descriptions of the sequences in SWISS-PR,OTfor the following strings: "SERISE/ persona1 communication) using the program PROFILESEARCH (Gribskov rt ul.. 1990). The THHEOSTSF:-PROTETSKTXASE. SER.;THRT'ROTEJS KTXXSR. PR.OTElS-S;ERTSE/ seven kinases used to generate the protile ;UT. 'I'HHEOSISI.: I<TS.WI<.PHOTk:TS-SER.,WHR bovine cA3rIP-rlependent protchin kinast. (I'IR wdc 1516 zyxwvutsrqpo zyxwvu Hidden Markov M.odels zyxwvutsrq zyxwvu zyxwvuts zyxwvuts zyxwvuts KINASE, . TYROSINE-PROTEIN KINASE, TYR-PROTEIN KINASE, PROTEIXTYROSINE KINASE, PROTEIN-TYR KINASE, V-ABL, C-ABL, V-FGR, C-FGR, V-FMS, C-FMS, V-FPS/FES, 1’-FES/F’PS,C-FPSFES, C-FESFPS, V-FYN, C-FYN,V-KIT.C-KIT, V-ROS,C-ROS, V-SEA, C-SEA, V-SRC,. C-SRC, V-YES, C-YES,, V-ERBB”. Of the 296 SWISS-PROT 22 sequences that were above the Z-score cutoff of 6.0 and were thus classified as containing a kinase domain b y our HMM, 278 were similarly classified by PROSITE, PROFILESEARCH and the keyword search. These 278 sequences may be considered to constitute “certain kinases”.Figure 11 shows the multiplesequence alignment generated by our HMM of some representative kinases from thisset (sequences 1 to 22). Sequences 23 to 40 are the 18 sequences (296 minus 278) that were designated as kinases by the H M M and one or two of the three othermethods. For PROSITE, we consider a sequence to be a kinase if it satisfies one or more of thethreepatterns PROTEIS-KINASE-ATP,PROTEIN-KINASEST or PR.OTEIN,KINASE-TY R as a true positive (“T”in Fig. 11B).PROSITE false negatives (“N”), potential hits (“P”) and false positives (“F”, Hidden Markov Models zyxwvu zy 1517 zyxwvutsrq zyxwvutsrq zyxwvut zyx zyx sequenceswhich do not belong tothesetunder consideration) are ignored. Among the 18 sequences classified as kinases by our HMM, eight (23 to 26, 35, 38 to 40) were also deemed to be kinases by the keyword search and PROSITE, andone (27) byPROFILESEARCH and PROSITE. The remainder (28 to 34, 36 to 37, those indicated by yo in Fig. 11B) are particulate guanylyl cyclases and except for 36 to 3 i , PROFILESEARCH also defines them as possessing a kinase domain. These guanylyl cyclases contain a single transmembrane domain, a cyclase catalytic domain andan i n t r a r ~ l l ~ ~ l protein ar kinase-like domain inwhich protein kinaseactivity has not beenseen to date (reviewedby Garbers, 1992). Although these sequencesare not kinases in terms of function, they possess all the conserved subdomains (subdomain I, the nucleotide binding'loop is modified in some) and the majorityof conserved residues present in certain kinases (see Subdomain of Fig. 11A and positions indicated by *). Sequences 41 to 50 are the top ten sequences in SWISS-PROT immediately below our cutoff of 60. Of these, the first three (41 to 43) were classified as kinases bytwoout of PROSITE, PROFILESE.4R.CH and the keyword search. Our cutoff was zyxwvutsrq zyx zyxwvu zyxwvu zyxwvuts zyxwvuts zyxwvuts chosen from a visual inspection of a histogram of Z-scores which indicated that 6 0 lay in a large gap (see Fig. IO). If the 2-score cutoff is lowered to the next largest gap (from 2-score 3.9 to 4-8)between sequences 43 to 44, then these t h r e e viral sequences (41 to 43) would also be categorized as kinases by the HMM. Of the eight sequences (41, 51 to 53, 56 to 57. 59 to 6 0 ) that were not classified as kinases by our HMM but wereclassified only by the keyword search and PROSITE, one (41) is t h e first sequence below OUT cutoff discussed above. Four (56 to 57, 59 to 60) are partial sequences where the kinase domain is ahsent. Three (51 t o 63)possess divergent forms of many of the cvnserred regions and like SI to 43. although they are below our cutoff, the HJIM is able to generate an alignment that correctly identifies divergent forms of conserved regions. FinaIly, there are three aminoglycoside 3'-phosphatransferase sequences (5.4 t o 55,58) which are only designated as kinases berause they sat.isfy the PROSTTE espression for the catalytic loop. Inspection 0f'Figur.e 1 1 B permits an estimation of the accuracy of the various methods in distinguishing kinases from non-kinases in database discrimination tests. The HMM generates six false Hidden Markoa Models zyxwvu zy 1519 zyxwvutsr zyxwvu zyxwvutsr zyxwvu negat.ives (41 to 43,51 to 53) of which the first three fall immediately below our kinase cutoff. For PROFILESEARCH,thereare 12 false negatives (23 to 26, 35, 38 t o 41, 51 t o 33) but it should be recslled that eight of these (those indicated by $ in Fig. 11B)do not appear in the results obtained from searching SWTSS-PROT 25 provided to us bx W. Gribskov (personalcommunication). 1Ve suspect t h a t at least four (23 to 26) would be correctly classified as kinases by PROFTLEFEARCH leaving an estimate of three to eight false negatives. In the case of PR.OSTTE, using 0111- assumption of it kinast3 to be a true positive ( T )s q w n w I i w m y I!W til' three patterns, there are three false negatives (39,42 to 43). However, the 8ctua.l performance of t h e PROSITE patternst.hemselvesis much worse; scans of SWISS-PROT 22 with each of thepatterns PROTEIXKIXASE-ATP? PROTEI'?;KISASLf5T and PROTEIKXISASE-TYR individually yield 4 0 , 2 and 3 false negatives, respectively. The difficulty in quantifying the precise number of false positives and false negatives produced by the databasediscrimination tests may beillustra.ted l y emlhyinp an dternutire mechanism for assvssil1c t h 11~1nkwrc)f false ntyativrs. if sin1pIy zyxwvu 1520 zyxwvutsrqpo zyxwv Hidden Markov Models zyxwvutsrqpon zyx zyxwvutsr zyxwv Figure 11. A, Multipte sequence alignment generated by our kinase HMM of some of the sequences used to train the HMM (1 to 22) and test sequences from the SWISS-PR.OT 22 database (23 to 60) (see Results section (b)). Pu'umerals appearing in the alignments indicate the number of amino acids to be inserted at that point. otherwise the notation follows the convention of Fig. 5. I n Subdomain, the Roman numerals and * refer to the subdomains and residues conserved across 75 serine/threonine kinases given by Hanks b- Quinn (1991). A and B in PROSITE refer to the ATP binding and catalyticregions, respectively, used to create 2 different signature patternsfor kinases. X-rayidentifies the location of the mhelices AA-AI and /3-strands Bl-B9 (read vertically) derived from the 2-7 d crystal structure of the 1) (Knighton d a/.. 1991). Sequences 1 to 22 are catalyticsubunit of cAMP-dependentproteinkinase(sequence representative kinases taken from the March 1992 Protein Kinase Catalytic Domain Database (Hanks & Quinn, 1991). These are: CAPIC-ALPHA, cAMPdependent protein kinase catalytic subunit. a-form: WEE1 +. reduced size at division TTK, mouse serine/threonine kinase: S P K l . S.cerekier kinase cloned with anti-pmutant wild-type allele gene product; Tyr antibodies; RSKl-E, amino domain of type 1 ribosomal protein 66 kinase; PYT,putative serine/threonine kinase cioned with anti-p-Tyr antibodies; PKC-ALPHA, protein kinase C, z-form; PDGFR.-B. platelet-derived growth factor receptor B type; PBS2, polymix in B antibiotic resistance gene product: MIKl. S. pombe mikl aets redundantly with weel +; MCK1,S.cereviaiae protein kinase; 1XS.R. insulin receptor: HSI'K. Herpes simplex virus-US3 gene product; ERK1. rat irisulin-stimulated protein kinase; EGFR, epidermal growth factor receptor {cellular homolog of v-erbB); ECK, receptor-like tyrosine-kinase detectedin epithelial cells; D P T K l . developmentally regulated tyrosine kinase in D. dkc.&ietLm; CLK. mouse serine/threonine/tyrosinekinase; CDC2HS. human functional homolog of yeast cdc2+/CDC28; of u-src; and CAMII-ALPHA. calcium/c&nodulindependent proteinkinase 11. 1-subunit: C-SRC.cellularhomolog C-RAF, cellular homolog of v-ruf/mil.Sequences 2 to 4.6. 10, 11, 14. 17 and 18 are t.he candidate dual-specificity protein kinases as defined by Lindberget al. (1992). Sequences 23 to 40 are the SWISS-PROT 22 sequences designatedas kinases by our HMM {Z-score >6.0) but not by all 3 other methods, PROSITE. PROFILESEARCH and the keyword search. Sequences 41 to 50 are the top 10 sequences below our cutoff of 6 0 and 41 to 13 and 51 to 60 are sequences that were not classified as kinases by the H M M but were so by one or more (but not all) of the 3 other methods. Kote that sequences identified as kinases by all 4 methods are not shown. All sequences that are less than 200 residues in length zyxwv zyxwvutsr zyxwvu zy zyxwvutsr zyxwvutsrq zyxwvutsr zyx Hidden Markov Models the number of sequences denoted as kinases only by all three other methods is evaluated, the number of false negatives for each of the techniques differ from .for .the HMM the more detailedanalysis:two (42 to 43), seven for PROFILESEARCH (23 to 26, 35,38,40)and none for PROSITE (ignoring known false negatives as above). This general problem is further highlighted by the guanylyl cyclases (indicated by yo in Fig. I 1 B}.If the definition of a kinase is based upon function and not possession of particularsequencepatterns,thenthe guanylyl cyclases are the only false positives for both the HMM andPROFILESEARCH.ThePR.OSITE patternsPROTEINXINASEATP, PROTEI?;KINASEAT PROTEIN-KISASE-TYR and produce eight, none and two false positires, respectively, giving some indication of the act.ual PROSITE performance. OveraIl, both the HMM and PROFILESEARCH appear to perform generally better than PROSITE in the discrimination tests, with the HNM possibly having a slight advantage over PR.OFILE- SEARCH. The HMM database search did not suggest any new putative kinases in SWISS-PR.OT 22. However, a comparative examination of the HMM producedmultiplesequencealignmentand the crystal structure of the catalytic subunit of CAMPdependentproteinkinase(Knighton et al., 1991) (sequence l), a template for the protein kinase family, yields insights into the consen-ed regions and their functionsin kinases of unknown structure. Figure 11A dispIays the location of secondary structure elements obtained from this crystal structure. An invariant Asp in subdomain., VIb (Asp166 in Knighton et d.,1991) that is proposed to be the catalyticbase is known to diverge in guanylyl cyclases (28 to 34,36to 37) even though the immediate region is highlyconserved (Garbers, 1992). Our results indicate that other invariant residues appear to bereplaced as well. In the sea urchin spermatozoan cell-surface receptor for the chemotactic peptide “resact” (sequences28 and a), a Lys 1521 in subdomain I1 (Lys72) that forms part of the ATP a- and p-phosphate binding site is changed to His. Theheat-stableentertoxinreceptor of rat (36) replaces an Asp in subdomain IX (Asp2OO) that contributes directly to’stabilizationof the catalytic loop by Glu. Yeast VPS15 (sequence35), a probable serinelthreonine kinase that is autophosphorylated, lacks many of the residues in subdomain I. I n addition, a conserved ion-pair that stabilizes ATP (GIu91-Lys72) would be disrupted in VPS15 because the Glu in subdomain 111 is altered to Arg resulting in the appositionof two positively charged residues. In the putativeB12 kinases of two strains of vaccinia virus (42 to 43),. the proposed Asp catalytic base is replaced by Lys (cf. guanylyl cyclases). This is accompanied by a further change in the “general“ sequence of the catalytic loop: the normally positively charged residue a t R 2 has been altered to Glu. In genera1, all the sequences below our cutoff and the lastone above it (40 to 60) appear tolack s-helix F ( s e e X-ray in Fig. 11A). The functional and or structural consequenoes of these modifications on any kinase activity are not clear. ’ + (c) EF-hand experiments zy For theseexperiments we used the June 1992 database of EF-hand sequencesmaintainedby Kretsinger and co-workers (Nakayama et al., 1992). Sequences in this database are proteins containid one or more copies of the EF-hand motif, a 29 residue st.ructure presentin cytosolic calcium-moduIated proteins (Nakayama et al., 1992; Persechini el al., 1989; Moncrief et al.,1990). These proteins bind the second messenger calcium and in their active form function as enzymes or regulate other enzymes and’ structural proteins. The motif consists of a n a-helix, a loop binding a C a 2 + followed by a second helix. Although a number of proteins possess the EF-hand motif, some of these regions have lost their calcium-binding property. For our training set, we extracted the EF-hand structures from each of the ‘242 sequences in the f have been removed. B, Details on sequences 23 to 60 shown in the alignment (arrangedin order of decreasing 2-score). PrZL-score snd 2-scoreare measures of how well the kinase HMM fits these SWISS-PROT22 test sequence that were not present in the training set [ s e e Results section (b)for more details). I n HMM,PROFILESEARCH and Keyword, + denotessequences that are classified as dontaining a kinase domain and - those that do not. For PROFILESEARCH, -$ identifies sequences that do not appear in the results obtained’from searching SWISS-PROT 25 (not 22 as in HMM. Keyword and PROSITE) provided to us by X . Gribskov (personal communication). TwoPROSITE signature patterns for eucaryotic protein kinases havebeen derived andthese are labeled A and B in the alignment. A is the region believed to be involved in ATP binding (PROSITE entry PROTEIKXIXASEATP) while I31 and B2 indicate the area important for catalytic activity in serine/threonine kinases (PROTEIN-KINASE-ST) and tyrosine kinases (PROTEIEi-KIXASE-TYR), tespectively. In all instances, T signifies a true positive; F a false negative (a sequence which belongsto the set under considerationbut which is not picked up by the pattern); P a “potential” hit (a sequence that belongs to the set but which is not picked up because the region that contains the pattern is not pet available in the data bank, Le. a partial sequence); and ? an unknown (a sequence which possibly could belong to the set). * Indicates SWISS-PROT files which contain a cross referenceto the specified PROSITE pattern, but these PROSITE entries do not contam a correspnding pointer to the SWISS-PROT file. - Signifies sequences t h a t do not satisfy the kinase pattern and yodenotes particulate forms of guanylyl cyclase receptorswhich contain an intntcellulsr protein kinase-like domain but which have not been shown to possess kinase activity to date (reviewed by Garbers. 1992). zyxwvuts zyxwvutsrq zyxwvut zyxwvuts database. obtaining 885 EF-hand motifs having an average length of 29. For our first. experiment we brained five HIIJfs on all 885 EF-hand motifs. using the standard techniques described earlier. (In subsequent experiments,described below, we trained on sma.ller subsets of these 885 sequences.) The best model had a final lengthof 29, and a NLL-score (the average -log P(sequenceImode1)) of 61.41. As described in Methods section (e), wemodified the final model to enable it. t o sea.rch the SWISS-PRO"database forsequences containing the EF-hand motif. We computed .Z-scores for all sequences as described in section (0 and Figure 12 shows the resultinghistogram. In contrast to the kinases,a visual inspection of the histogram of Z-scores did not indicate the presence of a distinct gapthus makingtheselection of a cutoff more difficult. After choosing by eye a cutoff of $75 and excluding d l sequences with unknown residues (Ss). the model classified 232 sequences as containing the EF-hand sequence motif. As with the kinase experiments in the previous section, false positives and falsenegatives were identified in the following manner. A list of "certain EF-hands" u-as created from the union of sequences determined t o be containing the EF-hand motif by threeindependent sources: PROSITE, a kep-ord search, and the results of Michael Cribskov's PROFILESEARCH. ,Details of the PROSTTE and keyword searches are given by Krogh et nl. (1993n). THO different PROF1LESEAR.CH experiments &ere conductedfor us by M. Gribskor (personal chmrnunication). The first employed a. profile generated using the multiple sequence alignment of sequences classified as EF-hands by our HNJI and the second was constructed using an alignment of the following four sequences: Escheriehin coli galactose bindingprotein [JGECG, 1 EF-handmotif). 2), human troponin rabbitpawalbumin(PVRB, (TPHUCS, 4) and human calmodulin (NC'HL'. -I). Although a sequence may possess multiple copies of t.he EF-hand (or any other) motif, only the one which most closely resembles t h a t described by the H M M is identified. Of the 232 SWISS-PROT 22 sequences that were above the cutoff (Z-score >475) and were thus classified as containingan EF-hand motif by our HMMM,163 were similarly classified by PROSITE, bothPROFILESEAR.CH experiments and the keyword search (if only one of the PROFILESEARCH experiments isconsidered. then there are an additional 14 sequences making a total of 177).These may be considered t o constitute cekain EF-hands and Figure 13 shows t h e multiple sequence alignment generated by ourHMM of some representative EF-hands from this set (sequences 1. to 27).Of the 69 (232 minus 163) or 55 (232 minus 177) sequences above the cutoff and not categorized as EF-hands by all thr&other methods, 33 possess the mot.if but do not bind calcium (indicated by oh in Fig. I3B) and six ( 6 4 , 72, 88, 89, 91. 94) were classed as EF-hands by only one other method. The identification of certain EF-hands as compared t o certainkinases is notas st.raight- S W r d h h t h Figure 12. Histogram showing the numberof sequences w i t h a certain Z-score relatire to thp EF-hand model. zyxwvutsrq zyxwvut zyxwvu forward, making it difficult to ascertain the precise number of classification errors made by each technique. This problemarisespartly because of t h e absence of a pronounced gap in the histogram of Z-scores and the result.snt uncert.ainty in assigning an esact Cutoff (Figs 10 and 12). The mnemonic developed to identifF EF-hand homologs and distinguish them from a.nalogs (Sakayama et al., 1992) is known t o generate errors and is unable t.0 detect 8 of the P i sequences known to be EF-hands (sequences 1 to 27 in Fig. 13).Therefore. the sensitivity and specificity of the EF-hand database discrimination tests is unlikely to be comparable to the kinases. C'sing Figure 13H,an estimate of the false negative rate for each method was determined by using the simple notion of evaluating the number of sequences classified as EF-hands by a.ll methods other than the one being considered. (Those which possess the motif but do not bind calcium? denoted hy nb in Fig. 13R. are not considered.) Using this criterion. the number of false negatit-es are 1 for the HMJI (101). 20 for PROFILESEARCH using four sequences (28.47, 56 to 57,59.67. 74 to 82,84 to 85, 92 to 93. 96). seven for PROFILESEARCH using our EF-hand alignment (28,57,i4,79 to 80. 92 t o 93). one for the keyword search (58) andtwofor PROSITE (60,70).A similar analysis of fa.lse positives produces six for the HMM (52, 71: 83, 86. 90, 94) nine for PROFILESEARCH using four sequences (9i,99,11 1 to 112, I21 to 122, 129: 132 to 133):eight for t h e k e p o r d (123, 126, 130-131; 134 to 137) and one for PROSTTE (120). I t should be noted however, that a search of SWISS-PROT 22 using the PROSITE pattern EF-HASD produces different results: three false negatives and 24 false positives (compared with 2 and 1 using the simple criterion). A total of 26 sequences were not desig- zyxwvuts Hidden Markor .Models zyxwvu 1523 1524 zyxwvutsrqpo zyxwv zyxwvutsrqponmlk zyxwvut Hidden Markov Models .. .. .......... .......... .. .......... .... ..... .... .......... C.IItW.t.L D....SD...~C.UCL....I.P.IlLYI11~lgti~ Uk.rgllS?I L.PPII.E.1 D....TD...CSC.SISL....A.DI.LCl~-.pll~ll~~ ~,n(.rlR1F..l.ItpFL.l.l..........D-...~...CDC.KlCL....O~.~rlMI1K-~t~..b.iS9 ay.rdarv5'1L..Y.tElAE.L.A D FI...KDc.EtlV D . ~ . t 4 1 V p P - c b @ d ~ l O 4 m.ikCI69E l . W E . D . l . . . . . K 1t...LW..KVOI ....E.1F.VMAAI-lrikdb.dM g ~ i * I o I L 1.rNI.c.1 D....I~...G~.~~L....I.~.~~~~~~,~~~I gMqkIo5K I.?M.C.l.......... D....*5...c=.mL E.Q.OllM-#qddlp5J d'p~r4JU I.~al.I.1..........O....It...tDt.ll~....I.~.M~-~~~d~ IU+y(p7 F.11115.l .V.......... D....?t U C . W m . . . . P . ~ . l D ~ S ~ - ~ ~ * ~ ~ ~ 7 ~ wrph.&A..F.WU.I.L D. ..SI...Mt.EIDF O.~.oNcI-untilS m'i.~5751 .. .. .. .... .... ... ~ ~ . . .......... .... O.U..IKT~-rld.UiISO ..........0....SI...MI.EIDF ....P.R.ClTLSCl-u~119 ..........D ....OD ...IPC.mI....l . a . ~ - ? R - d r ~ ~ A y J n ................................ D....RI..,CDC.KlS~....D.R.~1~-~~~~~~25 U&.qllYI..L.OUII.I.I ..........O....TDD,..C%.SM ....I.P.lUCAR-.plIWIU4O pmqlimi8lU..I.tXU.K.l ..........0 ....U...IDI.KtSF ....K.R.OlFUP-migrdd&4 Ibld1da.d ..I.?MFL.I.f .._....... D ....M...CtC.fl- ....-.-.---U-a ......... r..iic ~ C . . L Z I.V I . L . 1 . . . . . . . . . . D . . . . W . . . C I C . U C L . . . . ~ . ~ . l ~ ~ l - ~ , ~ , j ~ 7 rag~~IlllE..F.UIYA.E.I ..........0 ....?E...ATC.LlXY ....L . O P . ~ 1 U - s p p l g t g i D t mUiW699E..ll.aKflA.A.A ..........D ....'to...ClG.tlIS ....T . D L . m - - I I - y i l U i . k I i -p*.p1400E ..F.KANA.I.1....._.... D ....PE...AIC.LlKY ....L . O V . ~ l q p p l @ g i 4 4 Yt.IkbIIL..F.IUFL.E.l urpI..a5JA..F.PILIt.t.L mllplAlllIE..WOIFI.E.L ..........0....AD...SIC.IIDl.... zyxw zyxwvut .......... C ....ll...CC.lILV ....E.CI.VCl9UL-t....&d12 .........D....lI...SDC.PLOF....O.W.-----.......... ..-4. ILLS. ...D. W.lAILICK-C4%.ap.. .. IllmrdYE. .L. pLL1.L.F. . ........-....--. r ~ . l . ~ t..nI .~ I ~ . S ........... A D....LD...KlC.-lsK....a.a.~~~-~~g~i~~7 u.+..q4Ol..L.LUAl.A.r .. ~ U . r t ~ l l1.DUII.L.L. ' ....... .. .... rlltl~l511 L.LCGFl.C.C... D....VE..-IDe.SlSS P.Fl,.TIutAS-pLp9dC247 ~rU~b%E..L.UALS.L.Lqpllplqll~ U . . . M % t U F ....E.P.ULIUL-qm.aldl.49 IE...Lzc.-YSL ....1.1L.UIL1AA-gydlrk~JS m~ryup41L..I.NVA.I~~ l a ~ l U 5 I LK . I ~ . l . l . . . . . . . . . . I....VE...Ml.SlT P.OL.PsnolI-~h~qblS .... ..........- .... .. .. .... u ~ . W I ~ SL.LDA~.L.1..........D....OLir~C~.OIS~....O.l~.Sl~~-d~d~~~ 1llclql5lL L.LCCn.C.C D. IE...lffi.SVSS O.FL.lI~l-p~~kd~47 rrr.rtIOlY L.EL+n.A.L..........L....K1...~I~.~LDA....D.P.UItltL-~d.dcU116 .ir&i&ISl4E..F.KINA.E.Y D R...AIC.tIKR ....L.OV.rILLUI-qppl~@I6 .. .......... ... .. .... .......... .... ).-.D...I.EIY~I.F.EB~I~~ ....I.OV.rUlrsL-mkd.glo2 mingld21L..V.FtCFE.V.t .......... D....vI...IDc.C-DF Ugll).ld31L ..L.LAIL1.E.C B..~.~...IlC.U~S....O.Dl.~WU~-~~l~~~JOl r . l c l q l 5 I K ..L.UCFI.C.C...: D lE...IDC.5155 O.FL.lIUIlA-pl~p6~241 ~~Itlql5II..L.LCCCI.C.C D ....IE...lDC.SlSS O.FL.TILUlI-pkpk6t241 .rfpl.IO.S..Y.LIVL.P.T ..........D ....LD ...t S C . b J J ....V . f l . U A L F 5 A - ~ k l s n k l 6 7 0 ID...UC.TlEC K.FI.OKY-M-m*dqhvpGd -..-.-trYI.I.F ..I~.I~I~O~BII.TLT~L.L.C D r ~ . . m . u .... n ~L.EY.C~CFCR-W~~A..II* ~1.~~5~11-1Vrn.E.C 0.. GO. ..0 0 K . L l n ....I.W.CEffA~-..di"~lt .......... ...... .... .......... .......... .... .......... .... ......,.._ .... .... .... .......... .. .. g ~ U l . 2 ~ L 1.-~.U.~.~..........K....~...S~.KIU....-.--.---~-~~~~~~L~ -~.~~~..T.~~P.I.A..........D....-...~.R~....L.F;.~~~~~~~~~~ ..........IlrEC...CDA.WI....~lr.LD~-g7rkkql.l5 ..........D ....LO...tOK.VIIL ....F.W.ACQCI1-meImdl.i .......... - ....--...---. ? . S .... I l.O.LSIttIE-L~irk$16 .. ki.dfl(22SL L.CDfA.0.l ~~~iii2IIO~T.TVFE.t.C dgh 5..1.-----.-.rrrd~d.d.1)91 ..I . O I M . 1 . L arriiilBLbcT.lm.1.C ..._... ..........E . . . . S L . . . L D ~ I ~ C K . . . . C . E I . ~ ~ ~ - ~ q k ~ k ~ 2 ~ ..........D ..._LD...IDI.T:lL ....0.Fd.rccrCtK-qkdidkdlvi Plriii265Lb'T . ~l . C .......... . 0.. . . L D . i mi ........E..L.tEOlIiL.C ..........Dg ...lE...COL.CTU ....Y.~.llbUM-.gidtm69 zyxwvutsr ..........-....-...--. ..I . D A l l ~ . V . L L ~ . .S, . V . l L . ~ t I l E U ~ ~ l i ~ ~ ~ 5 0 .I~~~~~E..L.U~YD.L.L..........-. L . O T L ....L . K c . D O U U L - k i ~ ~ 4 3 y ~ l l l 5 CF.KC--C.t.l ..........s.... M~..COD.----....-.--.G5Il~-~~~i49 ~ r lfUa.. r l.M~.L.I.......... ~ l . . l D . . . C ~ . l I ~ . . . . l . I ~ , ~ l D - ~ - ~ p ~ ~ ~ Y.CC@cllD .. ..--...-- ~ ) *L..UYE.~.I..........C....tl...L~.CLCF....O.~.O1OA~-~~il~~qlIa .. ...__._... -....--...-s.nt1__.. ~.a.rarrmr-.~r~,.ul .gLLrpr5921..t.F.Ucc.1.s ..........s ....~ . a.c . .olc~'P.El.ct~~c~-l'~,.i'7~ I~ilqdl .wt..I.sIT..T.ms.5tyr*prl((O ..L.l~CS.R.Cbg~. .....-....~...~-. ~ : S ~ . . . . ~ . ~ . ~ ~ ~ - ~ = ~ b ~ ~ 7 6 A (cod! Fig. 13. nated as EF-hands by theHMM but were classified so by PROFILESEARCH, PROSITE or the keyword search. Of these, 19 were classified as such by only one of these methods. This includes five fragments where theEF-hand motifis missing: human and murine spectrin alpha- and beta-chains (123, 126, 131, 134) and rabbit calgizzarin (125). Inspection of the HMM produced alignment and examination of the putative calcium-binding ligands (Fig. 13) for the, 20 sequences immediately below the cutoff (97t o 1.16) and the false negatives and positives suggests that many possess potential EF-hand motifs. This includes six sequences whose 2-scores lie above our cutoff but are not classed as EF-handsbyanyother method: chicken myosin light chain alkali, smooth muscle (52); bovine calpactain I lightchain (71);Arabidopsis thaliana inorganic pyrophosphatase (83); rat placental calcium-bindingprotein (90) (note however that the mouse protein, sequence 88, is designated an EF-handbythe keyword search); and rat and bovine 1-phosphatidylinositol-iij-bisphosphate phosphodiesterase I11 (86 and 94). A notable example among the false negatives is the a-l subunit of L-type calcium channels from carpand rabbit skeletal muscle (97, 99) and rat and rabbit cardiac muscle (111, 112). These proteins play an important role in e.ucitstion-contrsction coupling and carry the calcium antagonist binding domains (reviewed by Grabner et at., 1991). They possess a highly conserved and evolutionarily preserved putativeintracellular region of I55 residues near the carboxyl terminus imrnediat.ety following the fourth internal repeat. This region has been suggested to ' Hidden Markor Models zyxwvu zy 1525 B zy zyxwvuts zy Fig. 13. . .. contain functional domains that are typical or essential for all L-type calcium channels regardless of whetherthe?; couple to ryanodinereceptors. conduct ions or hoth(Grahner cf nl.. 1991). The inferredEF-hands for these proteins owur within this cbonserred 155-residue segment. The above resultswere for an H3131 trained on all 885 EF.hnnd motifs ft-nln the Krt.tsjI1pr d n t a h a r . 1526 \ zyxwvutsrqpo zyxwv Hidden Markov Bodels zyxwvutsrqpo i .' B zyxwvu zyxw zy Figure 13. A, Multiple sequence alignment generated by our EF-hand HMM of some of the sequences used to train the HMM (1 to 27) and test sequences from the SWISS-PROT 22 database (28 to 137) (see Results section (c)). In Structure. H and L denote residues in an a-helical orloop conformation based upon EF-hands of known struct.ure(Xakayama rf d.. zyxwvu zy zyxwvutsrqp zyxwvutsr zyxwvutsr zyxwvuts zyxwvut Hidden Markov Models - There is considerable overlap between this training set and the EF-hand motifs found in SWISS-PROT 22, so in order to provide some clearer cross validation of ourresults we alsodidanother series of experiments.Intheseexperiments, modelswere estimated using training sets consisting of different numbers of randomlychosenEF-hand sequences from the database of 885 EF-hand sequences. For training sets consisting of 5 , 10, and 20 random EF-hand sequences, 15 models were estimated, each using a different randomly chosen training set. For training sets consisting of 40, 80, 100, 200, and 400 random EF-hand sequences, five models were estimated. In all, 70 models were estimated. A model's performance after training was gauged on how well i t performed on a test set which consisted of motifs from the database of 885 sequences that were not used in the training set. Thus for each model, two XLL-scores were computed (see Methods section ( f ) ) , one for the training set and one for the test set. These XLL-scores serve as a quantitative measure of hoG well the model is representing the sequence data. Figure 14 shows that for smalltrainingset sizes, the model overfits the training data. This is shown by low training NLL-scoresbutvery high testing XLL-scores. This effect largely disappears when the training set size reaches about 100 sequences. A model's performance was also gauged on how well i t searches a database for sequences containing the EF-hand motif. For each training set size, one model was randomly chosen to search SWISS-PROT 22. A histogram of theresultant Z-scores was plotted and a cutoff was chosen by eye. Thenumber of false positives was computed, as 1527 Tmining ssl Sam Figure 14. Average XLL scores for test and train sets for models with training sets of size 5 , 10: 20, 40, 8 0 , 100, 200 and 400.Error bars represent one standard deviation. described earlier in this section, by taking a list of certainEF-hands (i.e. determined to contain the EF-hand motif by the 3 independent sources) counting the number of sequences abovethe 2-score cutoff in the HMM database search that were not in thecertainEF-handlist.Figure 15 shows t h a t models built from smali training sets havelarge numbers of false positives. Again, this effect dis- 8.4 zyxwvutsrqponmlk 1992). PROSITE denotes the positions used to generate the pattern EF-HASD. Ca-binding identifies t h e 6 residues involved in octahedrally coordinating the calcium ion (denoted by X,5'. 2. x, z and y). The oxygen atom at. position y comes from the main-chain and so can be supplied by any amino acid. Sequences 1 t o 27 are representatives of the various EF-hand subgroups in the June 1992 database of EF-hand sequences maintainedby Kretsinger and co-workers (K'skaysma et al., 1992). These sequences are: CAMHS, Homo sapiens calmodulin; a A C . G , Callus gallus a-actinin; VISIKIX, G. gullua visinin; TPP24CF. Canis jamiliaria p24 thyroid protein; TPHUCS, H. sapiens skeletal troponin-C; TPAPl? l s t a c w ponlasliew troponin-C-1:TCBP25, Tetmhymena themophila TCBP-25; SPEMA, h70ngyloceutmha p p r d w speck; SCBPBLI. B r a n c h i o s W lanceolatum SARCl; QUIDLX?Loligo peulei squidulin; MOHSCR, H. sapiens myosin (RLC-ventricle); MOHSA1, H. sapiens myosin (ELC-LI-skeletal); LPSlA, Lylechirus p i d w a-Lpsl;. LAVl, Physamrn polycephulum LAW-2; EFH5, Tlypanosoma b m e i putative calcium bindingprotein; CVP, 3. Zance&um calcium vector protein; CRGHS. H. sapiens celmodulin-ielated gene; CMSE, Saccharopolysy eythraea bacterial-CAM; CDPK, Glycine tmxz calcium. dependent protein kinase: CDC31, Saccharomyces cererisiae ce11 division control protein 31; CALPLHS. H. sapiens calpain (light): CALCIB, Bos taurw calcineurin-B; CALBXGG, 0.gdZus caibindin; CALICE, Ctwmhbdab elegana cal 1 gene; BCHS, H.sapiens fi 5-100 protein; AEQAV1. rleqzwrea victuria aequorin-1; and 1F8,Try?lanosoma cruzi flagellar calcium binding protein. 28 to 96 are the SWISS-PR.OT 22 siquences designated as EF-hands by our HMM (2-score >475) but not by all 3 other methods, PROSITE, PROFILESEARCH and the keyword search. Biote that sequences identified as EF-hands by all 4 methods are not. shown. 97 to 116 are t.he top 20 sequences below our cutoff of 475; 117 to 137 are sequences that were not.classified as EF-handsby the HMM but were so by one or more (but not all) of the 3 other methods. B, Details on sequences18 to 137 shown in the alignment (arranged in orderof decreasing 2-score).XLL-score and 2-score are measures of how well the EF-hand H M M fits these database test sequence t h a t were not present in the training set (see Results section (c) for more details). Tn HMM. PROFILESEARCH. Keyword and PROSITE + and - denote sequences that are and are not. respctively. classified BS containing an EF-hand motif by the 4 specified methods. For PROFILESEARCH. Cribskov and HMM indicate results based upon profiles generatedfromfour EF-hand sequencesandour HMMM alignments. T. X. P and ? in PROSITE have the same meaningLIS in Fig. 11. ?& indicates sequences whichpossess an EF-hand motif bur do not bind calcium. might even be lowered further. However, there will be a limit on how small the number of available zy zyxwvutsrqp zyxwvutsr zyxwvuts zyxwvuts zyxwvutsrqpo zyxwvutsrqp NUlTJ2Wlr&inQni*rg5equarar Figure 15. EF-hand database search false positires for models trained with 5. 10. 20. 40. 80. 100. 200. 400 and 885 sequences. ~~ ~ ~ ~~ ~ ~ ~ ~~ ~~ ~~~ appearssubstantiall? when thetraining reaches a.bout 100 sequences. 'i, sequences can be if one hopes to obtaina reasonable model starting from a tabula rasa. We believe that the answer to the problem of small training sets is to add more prior knowledge into the training process. One way to do this is by starting with a better initial model. We have performed severalexperiments in which we have started with a model obtained from a small s t of aligned sequences, and thentrained the model further using a larger set of unaligned sequences. These will be reported in a future paper. We find that thistechniquecanoften give better results. This also suggests t h a t one application of HMMs may be in maintaining multiple alignments as the number of sequences in the alignment grows. Each time new sequences are added to a dataset of homologous sequences. we can begin with the HMM based on the alignment of the previous set of sequences, train it with the larger dataset thatincludes the new sequences. andthencreate a new multiple alignment for the larger dataset from this HMM. Xot only \vi11 the new sequences be included in the new alignment. but the alignment of the old sequences may be improved by utilizing the statistical information present in the larger dataset.t .\nother way to add more prior knowledge into the training process is to use a more sophisticated Bayesian prior. We are currently exploring the use of a prior on the probability distribution over the amino acidsin a match state of the model consistine of a misture of Dirichlet priors(8rown et ai., 1993). Vsing such a prior is like "soft-tying" the distributions in the states of the HJIM. By soft-tying we mean a combination of the idea of tying states (see e.g. Rahiner. 1989). inwhich thenumber of free parameters is reduced by having groups of states all sharing the same distribution on the output alphabet (the 20 amino arids in this case), and the idea of soft weight sharing from Sowlrin t Flinton (1992). i n which the regularizer (in this, case the prior for the tlistrilmtion of aminoacids) is also adaptirdy motlifird during learning. LVe have shown that this method can be used to estimate good KF-hand models using substantially fewer training sequenvez. Other types of more sophisticated priors can be obtained by switching from the alphal)st of the prirnar_v sequences to a different representationImed more on thestructural or chemical propertiesof the amino acids in the sequence. W e plan to explore these as wellIt is interesting to note that \{-e haveobtained quite gootl results i n multiple alignment and data. base searching without using any special weighting schemes to rnakr up for the statistical bias in OUI training sets (see e.g. Sibbald Br Argos, 1990), 01 employing Ihyhoffs matrix or any of its analog! (see e . g . IVatrrman. 1989) to take explicit rnutatior ~~ set size 4. Discussion X new method to model protein families using hidden Markor models has been introduced. The method is capable of tapping into the tremendous amount of sta.tistical information contained in many unaligned sequences from the Same family. For the cases of globins: kinases and EF-hands. t h e results have shown that by using this method. it is possible to obtain multiple alignments that mirror structural aIignments, having anlythe unaligned primary sequences as input. The results h a w also shown that the model can be used successfully in database searches for putative analogs of srquenc.es in a given protein family or domain. Finally. we believe that the model itself is a oaluable tool for representing the family or domain. The HYJI method we have proposed requires that many sequences be available from the family o r domain one wants to model. Since the nurnl~er of sequences in the protein databases i s growing rapidlj-, this may be less of a problem i n the future. but itwill always he a serious issue. Currently. only a relatively small number of sequences are avai1nt)le for most protein families and domains. For the globin family, we found that 400 sequenres is certainly sufficient. Preliminary., results indicate that 200 is enough, and even as few as 70 may sufice if they are chosen~'carefully from our database of 6-28 (70 chosen at random w i l l l~ nearly all a- and&chains).Ourexperiments using snlaller numbers of EF-handsequencesfor. training. as described in Results section ( 0 ) . show a similar zyx zyxwvutsr zyxwvu zyxwvu zyxwvut zyxwvu zyxwvut zyxwvuts zyxwvutsrqp probabilities between amino acids into account. 1 t of the PR.08ITE-indexed domains in a single long also remains to be seen whether or not incorporating protein, using the Vit.erbi algorithm. The remaining any of these extensions into the HMJl approach w i l l portions of the sequence could Ire marked as yield e\-en better results. “unknown“. While this would notconstitute a We also believe that some of the errors ma.tle I)? eompleteparse of the sequence. it would be very our HJIM models are due to the fact that these useful in providing some automatic annotation of models aresuboptimal, in the sense thatthrir new sequences, which is of critical importanceas the XLL-scores are not as low as they c~ouldbe. This is ra.te of growth of the protein databases continues t o because the E31 procedure is not guaranteed t o find accelerate. A related approach to proteinannothe globally optimal model for a given training set. tatiou is given by Stultz et al. (1993). and a related In other experiments, reported by Hanssler rt al. FlJIJI-based 1)SA parser for E. coli is described by (1993), we trainedan HMJI for.globins beginnink Krogh et ai. ( I 9936). with a nlodel derived from tlw Bashford et al. (1 ! M i ) A comparative examination of the HMM aJignment, and obtained a slightly lower SI.I.-scow produced kinasemultiplesequencealignmentand than an?- model from our esperilnrnts wing E:JI on thecqstalstructure or thecatalyticsubunit of unaligned training sequences (208 compiretl t o c.4MP-dependent proteinkinase(Knighton el al.. 210.3). Hence. w e know that 1511 is not locating the 1991) indicates a number of conserred residues in globally optimal model in this case. An important kinases of unknown structure that may be suitable open problem is to find a reliaMe way t o prevent for further experimental study (see R.esults section Ell from getting sturk and returning a suboptimal (I))). Results from our database discrimination tests solution. suggest the presence of an EF-hand calcium-binding Another issue is the adequacy of the hidden motif i n a highly conserved and evolutionary Markov model itself a s astatistical model of the preserved putative intracellular region of 155 sequence variation within a protein family. (‘learly residues in the 2-1 subunit of L-type calcium chanan H3IN provides at best n **firstortlrr” nmtlel of nels which playan important role in excitationsequence variation. There are many kinds of inter- cmtraction c*oupling(see Results section (c)). This actions in proteins that are not easily nodel led I)? region has been suggested to cont-ain the functional HJIJls. for esample. painvise correlations between domains that are typical or essential for all L-t.ype amino acid distributions in J)OSitionSthat are widely calcium channels regardless of whether they coup!e separated in the primary sequence. but close in the to ryanodine receptors, conduct ions or both. O h three-dimensional structure (see e.g. Klinger & EF-hand HMJI indicates the following proteins may Brutlag (1993)). It would he very raluable to hare also possess this motif: chicken myosin light chain more general models that incorporate such interalkali(smooth muscle), bovinecalpactain I light actions while still remaining computationally tracchain. -4rabidopeis trllaliunn inorganic pyrophosphatase. rat placental caleiam-binding protein ancl rat table. We are currently exploring the potential of one. model class of this type to capture the baseand bovine 1 -phosphatidylinositol-4,5-bisphosphate pairing in RSA fa.milies (Sakakibara ct ai.. 1993). phosphocliesterase 111. and hopeeventually to incorporate some of the Although there are many experiments left t o be features of these models into our protein models. done. based on our esperience. we believe that Finallv, we are encouraged by the quality o f the HJIJIs Hnd the E31 algorithmhavetremendous multiplesequencealignmentsgenerated by the potential in the area of statistical modeling of WoHMMs and the accuracy of the database seawhes. logicd macromolecules. Currently. most of this For example, the kinase H3lY is able to align potential remains to be rea.lized. correctly class TIT receptor tyrosine kinases which possess a domain that. differs from other receptor tyrosine kinasesby the insertion of a stretch Of 70 to We thank Peter Brown. Smen Brunak. Richard 100 residues (see the insertion between the 1)and E Durbin. Harry Soller. Martin Vingron. Don Momis. and helices in sequence 8, the /I-chain of the plateletMichael Zuker for valuable comments on this work. \-my derived growth factor receptor, in Fig. I ]A). k17ith special thanks to R,irhard Hughey for implementing our respect to thedatabase discrimination tests, we soft.ware on a MASPAR. parallel machine and to would eventually like to see HX’vIs built for all the YASPAR for providing computer time on their machine domains and families currently indexed by for some of these experiments.and very special thanks to PROSITE expressions. In many cases. HMMs for Michael Gribskor for running his PROFILFSEAR.CH subfamilies could be constructed automatically program for the kinases and EF-hands so that we muld using the method described in Methods section (d). rompare the results to those found with the HIIN. This work was supported by XSF grants CDA-9I16fW and Once this is done,this mightthen lead t o the TRT-9123692.OXSR grant SooO1+915-1162. XIH grant construction of a simple “protein language parser” number 03117129. and a grant fmm the Danish Satural using HJi3I.s. This parser could be constructed by Srin1c.r Research (‘ouncil. The full a1ignment.sand Z-xc~ore connecting all these individual H M h in parallel tahles dewribed i n this paper are arailabte in elertronic into a single large HMM with a global BEGIX and form. ancl (*allbe obtained by anonymousftp fronl Ex’1)state, and a transition from the ESD state ftp.rse.ucsc.edu. Our HMLI builclinp program and other back to the BECTS state. In principle. this parser tools (written i n C‘) will also be ma& availahlr frcm thr should be cqal)te of finding all n(vurrt*tl(*es of eacsh salne f t site. ~ ~ zyxwvu zyxwvutsrqpo zyxwv zyxwvut zyxwvu zyxwvuts zyxwvutsrqponmlkjih zyxwvutsrq Hidden Markov Models I530 References Xbe. S. 8- Warmuth. SI. (1990). On the compubational by complexity approximating of distributions probabilistic automata.In Proceedings of the 3rd Workshop O R C o m p u l a l i d Learning T h e o y . pp. 5'2-66. Morgan Kaufmann, Rochester. ST. Allison. L.. Wailace. C. S. & Yee, C. X. (1992). Finite-state models in alignment the of macromolecules. J . Md. E v d . 35, 77-89. Asai. K..Hayamizu. S. & Onizuka, K.(1993). HMM with In Proceedings of the proteinstructuregrammar. Hawaii Intermtiom1 Conference on System Scielrces. pp. 783-i91. TEEE ComputerSociety Press. Los Alamitos. C.1. Bairoch. X. (1992). Prosite: a dictionaryof sites and patterns in proteins. Nucl. Acids Res. 20.1013-2018. Baldi. P. 8: Chauvin. 1'. (1993). A smoothlearning algorithm for hidden Markov models. ,Veural Computation, in the press. Baldi, P.. Chaurin, Y.,Hunkapiller, T. 8 " h e . 11. X. (i993). Hidden Markov models in molecular biology: new algorithms and applications. I n Adzurnces in SeuraE Information Processing Systems 5 (Hanson. C o r a n 8- Giles. eds), pp. 7-47-72. Morgan Kauffmann Publishers. San Mateo?CA. Barton. G. J. (1990). Protein multiple sequence alignment Method8 Enzymol. and flexible patternmatching. 183, M 3 4 2 8 . Barton, C. J. 8: Sternberg. M. J. (1990). Flexible protein sequence patterns: a sensitive method to detect. weakstructuralsimilarities. J . X o l . Biol. 212 (2): 389-402. Bashford. D.. Chothia, C. h Lesk, X. 11. (1987). \ Determinants of a protein fold: unique features of the i. globinaminosequence. J . Mol. Bid. 1%: 199-216. Berger, J. (1985).Stalistical Decision Theory and Bayesian Analysis. Springer-Verlag, Xew Tork. Bowie. J . C.. Liithy, R.8 Eisenberg, D. (1991). A method to identify protein sequences that fold into a known three-dimensional structure. Science. 253: 164-170. Brown. X. P.. Hughey. R., Krogh. A.. Mian. I. S.. Sjolander. K. & Haussler, D.(1993). C'sing Dirichlet mixture priors to derive hidden Markov models for protein families. In Proc. Firsf Int. PonJ on tnklligenf Systemfor ;Molecular Biology (Hunter. L.. Searls. D.. Sharik. J.. eds), pp. 17-55. AhAI Press. Wash. D.C'. Cardon, L. R. & Storrno. C . D. (1992). Expectation maximization algorithm for identifying proteinbinding sites with rariable lengths from unaligned DX.4 fragments. J . Y d . Biol. 223. 159-110. Churchill, G. A. (l-989). Stochastic models for heterogeneous DSA sequences. Bull Xath Biol. 51. 79-94. Dempster. A. P.. Laird, X. M. 8 Rubin. D. 8. (1977). Maximum likelihood fromincomplete data ria the EM algorithm. J . Roy. Statisf. Sa. 5. 39. 1-38. Dickerson. R..E. & Geis. I. (1983).Hemoglobin: Stttuetrre. Function. EtaIuiion a d Pathology. Benjamin/ Cummings Pub. Co, Menlo Park, CA. Duda. R. 0. & Hart. P. -E.(1973). Ptzflem Cla.ssijicution and Scene Anulysis. Wiley. Xew York. Everitt, B. S. 8; Hand. D: J. (1981). FiniteMirture Dislribdions. Chapman and Hall, London. Feng, D. F. 8: Doolittle, R. F. (1987). Prc -si 'c sequence alignment as a prerequisite to correct phytogenetic trees. J . Mol. Pvol. 25, 351-360. Garbem. D. L. (1992). Guanylyl cyclasereceptors and their endocrine. paracrine and autocrineligands. cell, 71. 1-4. Chnan. S.. Bienenstock. E. & Doursat, R. (1992). Neural networksandthebiastvariance dilemma. N e u d Compdation. 4. 1-58. Grabner. SI.. Friedrich, K., Knaus. H.-G., Striessnig, J., Scheffauer. F.. Staudinger, R.. Koch. W. J., Schwartz. A. &. Glossmann, H. (1991). Cahium channels from Cyprinus carpi0 skeletal muscle. Proc. .Val. dead. Sci.. U.S.A.88. 727-731. Gribskov. 31.. Liithy, R. & Eisenberg. D.(1990). Profile analysis. Methala EnzymZ. 183, 146159. Hanks. S. K. & Quinn, A. M. (1991). Proteinkinase catalytic domain sequence database:identification of conserved features of primary structure and classification offamilymembers. Yethod8 Enzymol. 200.38-62. Hanks. S. K., Quinn. A. M. 8 Hunter. T. (1988). The protein kinase family: conserved features and deduced phylogeny of the catalytic domain.Scienee. 241. q2-5'2. Haussler. D. 8: Krogh. A. (1992). Protein alignment and at the conference SeuraI clustering. Presented Setworks for Computing. Haussler. D.. Krogh. A., Mian. I. S. & Sjolander, K. (1992). Protein modeling using hidden Markov Technical Report models: analysis of globins. N"CC-CRL-92-23 Universit? of California at Santa CA Cruz. ComputerScience Dept.. SantaCruz, 95064. Haussler. D.. Kmgh. A.. Mian, I. S. & Sjiilander, K. (1993). Protein modeling using hidden Markov models: analysis of globins. In Proecedings of the Hawaii Inlernutional Conference on System Sciences, IEEE Computer Society Press. Los Alamitos, CA. Hunter. T. (1991). Protein kinase classification. Methode E R Z I J ~ M200. O ~ .3-37. durka. J . 8: Slilosavljeric, A. (1991). Reconstruction and J . .\id. Ecol. 32. 'analysis of human Alugenes. 105-121. Klinger. T.B: Brutlag. D. (1993). Detection of correlations in tRSA sequences with structural implications. In First InternalimaI Conference on Inteiligenl Systems for .Ifdecular Biology (Hunter. L..Searls, D. B: Sharlik. d.. eds). A A A I Press. Menlo Park. Knighton. D. R.. Zheng. J.. Eyck. L. F. T.. Ashford. V. .I..Strong. S.-H.. Taylor. S. 8.d: Sowadski. J. M. (1991). Crystal st.ruc?ture of the catalytic subunit of c:\clit adenosinemonophosphate-dependentprotein kinase. Science. 253. 407414. Krogh. A . . Brown. SI.. Mian. I. S.. Sjolander. K. 8. Haussler. D. (199b). Hidden Markov models in conlputational biolokv: applications to protein CCSC-CRL-93-32 modeling. Technical Report I'nirersity ofCalifornia at Santa Cruz. Computer Science Dept., Santa Cruz. CA 95064. Krogh. -4..Mian, I. S. & Haussler, D. (1993b). A hidden E . coli DXA. Markovmodel that finds genesin TechnicalReport UCSC-CRL-93-33 Cnivemity of California at Santa Cruz.Computer Science Dept.. Santa Cruz. CA 95064. lander. E. S. & Green. P. (1987). Construction of multilocusgeneticlinkage maps inhumans. Proc. S a t . h a d . Sci., U.S.A. 8 4 , 23634367. Lawrencp. C . E. 8; Reilly. A . A . (1990). An expctation maximization (EM]algorithm for the identification anticharacterization of cammonsites in unaligned biopcl. m e r sequences. Proteins. 7. 41-5 I . Lindherg. R. X.. Quinn. A. 11. 8: Hunter. T. !1992) zyxwvutsr , zyxwvu zyxwvu zyxwvutsrqpo zyxwvutsrqp zyxwvutsrqp zyxwvuts zyxwvut 1331 Hidden Markov Jlodels Jhal-specificityprotein kinases: will any hydroxyl do? Trends Biochem. Sci. 17. 114-1 19. Liithy. R.. McLachlan. A. D. 8: Eisenberg. D. (1991). Secondary structure-based profiles: use of structureconserving scoring table in searching protein for structural similarities. sequence databases Proteins: Slnrcl. Funcl. Gene!. 10, 229-239. Moncrief. S. D., Kretsinger. R. H. 8; Goodman. M. (1990). Evolution of EF-hand calcium-modulated proteins. T.. Relationshipsbasedonaminoacid sequences. J . Mol. E v o ~30,5’22-562. . Sakayama. S., Moncrief. X. D. ?L Kretsinger, R. H. (1992). Evolution of EF-handcalcium-modulated proteins. 11. Domainsofseveralsubfamilieshave diverseevolutionary histories. J . 2 0 ! . Evot. 34, ’ 416448 S o d a n . S. (1990). Maximum likelihood competitive learning. I n Advances in Neural Infwmalion Prmessing Sysfems (Touretsky. D.. ed). vol. 2. pp. 574-582. Morgan Kaufmann. San Jlat,eo. CA. Sowlan. S. J . & Hinton, G . E. (1992). Softweightsharing. In Advances in Seural Information Processing System 4 (Moody. Hanson 8: Lippmann. eds). 3iorgan KaufTmann Publishers. San Mateo. CA. Persechini. A.. Moncrief, X. D. & Kretsinger. R. H. (1989). The EF-hand family of calcium-modulated proteins. Trends Weurosei. 12 (1 1), 462-467. Rahiner, L. R. (1989). A tutorial onhiddenMarkov models and selected applicat.ions in speech recognition. P r m IEEE. 77 (2). 25-286. Sakakibara. Y..Brown. Y.,Underwood, R.. Man. I. S. & Haussler. D. (1993). Stochastic context-free grammars for modeling RIGA. Technical Report UCSC-CRL-93-16 University of California at Santa Cruz. ComputerScienceDept., Santa Cruz, CA 95064. Sibbald, P. & Argos, P. (1,990).Weighting aligned protein or nucleic acidsequences t o correct for unequal representation. J . Mol. Biol. 216, 813-818. Stultz, C. M., White. J. V. & Smith. T. F. (1993). Structuralanalysis basedon state-space modeling. Protein Sei. 2, 305-315. Subbiah. S. & Harrison, S. C. (1989). A methodfor multiple sequence alignment with gaps. J . Mol. Biol. 209.534-548. Tanaka, H.?Ishikawa, M., Asai, K. & Konagaya, A. (1993). Hidden Markov models and iterativealigners. In First I n l e d i d , Conference on Indelligenl Systems for Ydecular Biology, . U A I Press, Menlo Park. Taylor. IV. R. (1986). The classification of aminoacid Conservation. J . Theoret. Biol. 119. 205-218. Vingron. 11. & Argos. P. (1991). Motif recognition and alignment. for many sequencesby comparison of dotmatrices. J . Mol. Bwl. 218, 33-43. Waterman, 11. S. (1989). Sequence alignments. In Nafhematieal Yeflrods for DXA Sequences (Waterman. 31. S., ed.). CRC Press, Boca Raton, FL. Waterman. Y. S. & Perlsitz, 1. D. (1986). Line geometries forsequence ‘comparisons. Bull. &fa& Biol. 46. 567-577. White. J. Y.,Stultz, C. M. & Smith, T. F. (1991). Protein classification by nonlinear optimalfiltering of aminoacid sequences. Unpublished manuscript. zyxwv ’\ Edited by F . Cohen (Received 4 January 1993; accepted 23 September 1993)