Hidden Markov Models

1. Hidden Markov Models As used to summarize multiple sequence alignments, and score new sequences

2. Time flies like an arrow, but fruit flies like a banana • The Markov Property: the state of a system depends only on its previous state. • Works very well for speech recognition: a sound is modified by the sound immediately preceding it. • For DNA and protein sequences, there is no reason to think that the subunit just upstream is the only influence: probably the whole neighborhood is important. • So, HMMs are just an approximation, but they work very well. • In a first order Markov model, the current event only depends on the immediately preceding event. Second order models use the two preceding events, third order models use the three preceding, etc. • as far as I know, all the HMMs for protein sequences that are readily available are first order models.

3. Hidden Markov Models • Position-specific substitution matrices have an important weakness: they don’t deal very well with gaps caused by small insertions and deletions. • A we know, gaps are part of most sequence alignments. • The standard solution these days: profile hidden Markov models (HMM). • Like PSSM, they are a statistical model of a multiple sequence alignment, but they easily incorporate gaps. • HMMs can be used to evaluate a query sequence: to determine how likely the query is to fit the model.

4. Hidden Markov Models • The basic problem is: given a series of observations, determine the most probable set of hidden states that generate them. • For matching a query sequence to a model, each residue in the query is assigned one of 3 states: it is a match/mismatch to a residue in the model, or it is deleted, or it is inserted • but, to understand HMMs, we need a simpler and more concrete example, taken from Wikipedia, which we can call Life in St. Paul. • I have a friend (Carol) in St. Paul who has only 3 interests in life: going for a walk, shopping, and cleaning her apartment. • St. Paul weather can be either sunny or rainy. Carol decides which of her 3 activities she will do based on the weather. • Every day she e-mails me to tell me what she is doing. (exciting!) • From the sequence of her activities, I infer the most likely weather for each day.

5. Terminology • The HMM for Life in St. Paul has 2 states: rainy and sunny. • These states are hidden: we don’t know them and can only infer them from Carol’s behavior. • For each set, there is a set of transition probabilities: what will happen the next day. • Which depends only on what happened today (the Markov property) • prob R -> R (stays rainy) = 0.7 • prob R -> S (goes from rain to sun) = 0.3 • prob S-> S (stays sunny) = 0.6 • prob S -> R (goes from sun to rain) = 0.4 • For each state, there is a set of emission probabilities: what Carol’s behavior will be under different weather conditions. • These behaviors are our observations: she tells us what she is doing each day. • if state is Rainy: • walk = 0.1 • clean = 0.5 • shop = 0.4 • if state is Sunny: • walk = 0.6 • clean = 0.1 • shop = 0.3 • A path is a sequence of states through the model: what the state is for each day. We are trying to determine the optimal (most likely) path that fits our observations.

6. State Diagrams • One way to draw state diagram is to have arrows going back and forth between states. Note that there is also a transition back to the same state. • This represents eternity: each day is in one of the other state.

7. State Diagrams • Another way is to move forward in time, with a new pair of Rainy and Sunny • states for each new day. I like this because it is a closer approximation to • progressing through a sequence, trying to get a match between the query • sequence and the model.

8. Basic Questions • We are interested in 3 basic questions: • Given the model and its parameters (the transition and emission probabilities), what is the most likely path that explains our observations? • This is done using the Viterbi algorithm. • What is the overall probability that the model fits our observations? • This is the most common use of HMMs in bioinformatics: determining the probability that the query sequence fits the model • Uses the forward algorithm, which is very similar to Viterbi • Given a set of observations, how can we accurately estimate parameters for the model? • This uses the Baum-Welch algorithm.

9. Viterbi Algorithm • Used to determine most probable path (series of states) given a set of observed data. • the forward algorithm is closely related and can be performed at the same time: it generates the overall probability of the observed sequence. • based on the idea that if you always take the most likely next step, you will get the most likely overall path. Given a first order model, this seems quite reasonable. • Very similar to dynamic programming • Go through the observations from beginning to end, inferring a state of the hidden machine for each observation. • Keep track of 3 values: overall probability (prob), Viterbi path (v_path), and Viterbi probability (v_prob) (i.e. the probability of the observations given the Viterbi path). • Test all possible next steps. The probability of a possible step given its corresponding observation is the probability of the transition to that next state times the emission probability of the observation in that state. • for overall probability prob, multiply each new prob by the old prob, then add together (i.e. overall probability is sum of all possible next steps) • for Viterbi prob, just take the highest next-step probability and multiply it by the previous v_prob • for Viterbi path, add next state to list of previous states

10. Our Example • observations: WCSW • start at beginning. We will simply assume that the day before the beginning was sunny Observation is W (walk), so go with the most probable. • P(Walk | Sunny) = 0.6; P(Walk | Rainy) = 0.1 • prob (overall probability) = 0.6 • v_prob (Viterbi probability) = 0.6 • v_path = S • Step 2 (next day) observation is C (clean) • Given that previous day was sunny: • stays sunny. Transition prob = 0.6. Prob of cleaning on a sunny day (emission prob) = 0.1. Total probability for this possibility (i.e. it stays sunny) = 0.6 x 0.1 = 0.06. • turns rainy. Transition prob = 0.4. Prob of cleaning on a rainy day = 0.5. Total prob for this step = 0.4 x 0.5 = 0.20. Therefore this is the most likely state of step 2 • update overall probability: prob = 0.6 x 0.06 + 0.6 x 0.2 = 0.152 • update Viterbi probabilty: v_prob = 0.6 x 0.2 = 0.12 • update Viterbi path: v_path = SR

11. More • step 3. Obs = S (shop) • given previous day was Rainy: • stays rainy = 0.7. Prob of shopping on a rainy day = 0.4. Prob for this step = 0.7 x 0.4 = 0.28 • turns sunny = 0.3. Prob of shopping on a sunny day = 0.3. Total = 0.3 x 0.3 = 0.09 • best is R. v_path now is SRR • overall prob = previous prob (0.152) x sum of this step’s possibilities: 0.152 x (0.28 + 0.09) = 0.0562 • v_prob = previous v_prob x best current prob = 0.12 x 0.28 = 0.0336.

12. Still More • Step 4 (last). Obs = W • previous day was R: • R-> R transition prob = 0.7. Emission probability of a walk on a rainy day is 0.1. Total for this step = 0.7 x 0.1 = 0.07. • R -> S prob = 0.3. Emission prob = 0.6. Total = 0.3 x 0.6 = 0.18. • Best is S (turns sunny). v_path = SRRS • overall prob = previous (0.0562) x sum of current possibilities (0.07 + 0.18 = 0.25) = 0.562 x 0.25 = 0.0141. • v_prob = previous (0.0336) x best (0.18) = 0.00605. • We have now determined the best set of states that would give us the observed data (v_path = SRRS), and the overall probability of the observed data (given the model) (prob = 0.0141). • It is obvious that multiplying probabilities soon gets to very small numbers, so in general they are converted to logarithms and added.

13. Results from a Longer Run activity: CCSCCCSCWCWSCSCSSSCSSSWWWWSCCSSWCCSWWSCWSCSSCSSSWW weather: RRRRRRSSSRSSRSRSRRRRRRSSSSSSRRRSRRRSSSRSSRRRSSRRSS predicted: RRRRRRRRSRSSRRRRRRRRRRSSSSSRRRRSRRRSSSRSSRRRRRRRSS • Viterbi prob = 5.48e-30 • Overall prob = 1.92e-24 • Percent correct = 86.0 %

14. Estimating Model Parameters • The parameters that need to be estimated are the transition and emission probabilities. Since they must sum to 1, our St. Paul weather HMM needs only 2 transition probabilities and 2 emission probabilities estimated. • Sequence models have very large numbers of parameters to estimate, which requires large amounts of data. • If you have training data with known states (e.g. you know when it is sunny or rainy, plus your friend’s activities, for a period of time: parameters come from simply calculating the proportion of each transition or emission. • Example: on rainy days, Carol took 6 walks, cleaned 30 times, and shopped 24 times. Emission probability estimates are 6/60 = 0.1, 30/60 =0.5, 24/60 = 0.4. (very frequentist!) • problem of 0 counts. Add a “pseudocount” so this probabilty isn’t zero (because it might actually happen). Typically add 1 to numerator and total number of possibilites to denominator. • e.g. for walks, emission prob = (6 +1) / (60 + 3) = 7/63 = 0.111

15. More Parameter Estimation • If you don’t know the underlying states: you just have a series of observations: you can use the Baum-Welch algorithm to estimate parameters. • Need the backward algorithm in addition to the forward algorithm. The backward algorithm works exactly the same as the forward algorithm, except that you calculate from the end to the beginning. • The Baum-Welch algorithm starts with a set of parameters that have been estimated or guessed. • Then feed each training sequence through the model, with overall probabilities calculated both forward and backward. • Then use these probabilities to weight the contribution of each sequence to the transition and emission probabiliites. • Then use this to adjust the parameters.

16. Profile HMMs • Start with a multiple sequence alignment: we want to make a statistical model of it.

18. Building the Model: Match States • If we were just trying to do an ungapped alignment, we could use a very simple, unbranched HMM, where every match state can only transition to another match state • Each match state has an emission probability for each amino acid, with different emission probabilities for every match state • This is essentially what a PSSM is: the scores in each column of a PSSM could be scaled to 0-1, giving emission probabilities. • The transition probabilities are all 1: the only choice is to move to the next match state.

19. Insertion States • In a multiple alignment, there are often columns which are gaps in most sequences but have amino acids in a few sequences. • These columns are best dealt with as insertion states. • As you move through the model, emitting the query sequence, insertion states emit the extra amino acids found in these columns. • Insertion states have emission probabilities that are usually the same as the overall proportion of each amino acid in the database. • Insertion states loop back on themselves, meaning that multiple residues can be inserted here. • Insertion states are entered from one match state and then exited to the next match state: insertion occurs between adjacent amino acids.

20. Deletion States • Deletions are places in multiple alignments where most of the sequences have amino acids but a few of them have gaps. • Deletion states are used to skip over match states. • You can skip multiple match states by moving from one deletion state to the next. • Deletion states act like affine gap penalties: the transition probability for entering a deletion state from a match state is equivalent to the gap opening penalty, and the transition probability from one deletion state to the next is the gap extension penalty. • Unlike match states and insertion states, deletion states are silent: they emit nothing.

21. The Full HMM • A full HMM model for a making a profile of a multiple alignment. • Note that it is necessary to have begin and end states too, and often a bit more complicated than is shown here. These are needed to deal with local alignments as opposed to global. • HMMs also usually have transitions directly from insertion states to deletion states, but these are considered very improbable, but their existence helps in model building.

22. Determining the Parameters • Start with a good multiple alignment, which generally involves hand refinement. This is often refered to as “curated”. • Most of the columns will become match states. • Pseudocounts are needed to make sure that each match state has a non-zero probability of emitting all 20 amino acids. • Insertion and deletion states are given some background low transition probabilities to start with. • columns that are mostly gaps are considered insertion states, and gaps in columns with mostly amino acids are modeled with deletion states.

23. Using Profile HMMs • The query sequence is used as a set of observations that are emitted from the HMM if the proper path is followed. • Either the Viterbi algorithm (to give the probability for the most likely path) or the forward algorithm (to give the probability of emission by all possible paths) can be used. • To get log-odds scores, the probability of emitting the sequence from the model is compared to the probability of generating it randomly given the background frequencies of each amino acid.

24. HMM Collections • The J. Craig Venter Institute (formerly known as TIGR) has a large collection of HMMs, called TIGRFAMS. • They are classified according to their isology: how specific they are. The least specific are motifs, then domains, then families (super and sub), and most specific are equivalogs. • Equivalogs are homologous proteins which all have the same function. They are usually full length proteins, not just domains. • A good hit on an equivalog is strong evidence that you have determined the protein’s function. • Lots of useful information about the HMMs also. • Two cutoff scores: above the trusted cutoff there should be no false positives. Below the noise cutoff the match is of no value. • http://www.jcvi.org/cms/research/projects/tigrfams/overview/ • A bit slow, in my opinion. • Search results are not linked: you need to do the sequence search, then go to the text search to see what they mean.

25. PFAM • PFAM: probably stands for “Protein Families”. • http://pfam.sanger.ac.uk/ • Two types: PFAM-A is manually curated, and PFAM-B is automatically generated from a larger database. • PFAM HMMs are generally domains as opposed to full-length proteins. • The JCVI site also includes PFAM HMMs, but the actual PFAM site has search results with links to information about them.