EM with Many Random Variables Another Example of EM Sequence Alignment via HMM Lecture # 10

1. EM with Many Random VariablesAnother Example of EMSequence Alignment via HMM Lecture #10 Background Readings: chapters 11.6, 3.4, 3.5, 4, in the Durbin et al., 2001, Chapter 3.4 Setubal et al., 1997 This class has been edited from Nir Friedman’s lecture. Changes made by Dan Geiger, then Shlomo Moran. .

2. EM for processes with many dice In the previous class we presented the EM algorithm for the case where the parameters are probabilities associated with a single “die” (i.e., probability space/random variable) . In practical applications the model may include many dice (like the HMM model). The generalization of the EM algorithm to many dice is rather straightforward, and given next.

3. EM for processes with many dice The model is defined by the parameters (random variables, or dice) and the simple events. Let the random variables be Zl (l =1,...,r), Zl has mlvalues zl,1,...zl,mlwith probabilities {qlk|k=1,...,ml}. Each simple event y corresponds toa sequence of outcomes (zl1,k1,...,zln,kn) of the random variables used in y. Let Nlk(y) = #(zlk appearsin y).

4. EM for processes with many dices Similarly to the single die case, we have: Define Nlk as the expected value of Nlk(y), given x and θ: Nlk=E(Nlk|x,θ) = ∑y p(y|x,θ) Nlk(y), Then we have:

5. L(λ) for processes with many dices Nlk

6. EM algorithm for processes with many dice Similarly to the one dice case we get: Expectation step Set Nlkto E(Nlk(y)|x,θ), ie: Nlk= ∑y p(y|x,θ) Nlk(y) Maximization step Set λlk=Nlk / (∑k’ Nlk’)

7. EM algorithm for n independent observations x1,…, xn : Expectation step It can be shown that, if the xjare independent, then:

8. EM – One More Example The DNA of species in planet Melmek contains two letters: A and B. In the evolutionary process in Melmek, a mutation of a letter occurs in two steps: first the letter may be deleted, and in case the letter is not deleted it can be changed to the 2nd letter. The unknown probabilities of these events are respectively There are two species in Melmek, S (for Son) and its direct ancestor F (Father). Scientists were able to deduce that the DNA of F contained a sequence of two letters "AX", where "X" is equally likely to be A or B (Prob(X=A) = Prob(X=B) = 0.5). • Describe the probability space defined by the evolution of the two letters "AX" in F into a sequence (of up to two letters) in S: • Write down all the "simple events". • For each event write its probability (as a function of ), • and its contribution to the six statistics used by the EM algorithm.

9. Example (cont) • For writing the simple events, we assume the following order of “dice tossing” : • decide if the initial sequence is AA or AB (with probability 0.5 each) • Delete/replace the right letter (one or two tossing) • Same for left letter. • For instance, a simple event of deletion of two letters is:

10. Example (cont) • Two more simple events which start with AA: There are altogether 18 simple events: 9 which start with AA and 9 which start with AB

11. Example (cont) Later information showed that the sequence "AX" of F evolved into the sequence "A" in S. For given , write down the probability of the above scenario of evolution (of "AX" in F evolving into "A" in S).

12. Example (cont)  3. Write a single round of the EM algorithm to estimate the parameters which maximize the likelihood of the above scenario, starting with arbitrary initial parameters Show that the outcome is independent on the initial parameters. Calculate the counts of each outcome in each simple event: Regardless the initial parameters: del and remainhave exactly the same count, hence they receive each probability 0.5 Also, A→A and B→A both get probability 1.

13. Example (end) What are the parameters which maximize the likelihood of the above scenario? Justify your answer. Solution 2: The above is parameters are obtained after each iteration of the EM algorithm. By the EM correctness theorem, this must be the (unique) maximum.

14. EM for other discrete stochastic processes The EM algorithm is applicable to a general scenario in which we wish to maximize p(x|)=∑yp(x,y|). Where the experiment (x,y) is generated by a general “stochastic process”. The only assumption we make is that the outcome of each experimentconsists of a (finite) sequence of samplings of r discrete random variables (dices) Z1,...,Zr , each of the Zi ‘s can be sampled few times. This can be realized by a probabilistic acyclic state machine, where at each state some Zi is sampled, and the next state is determined by the outcome – until a final state is reached.

15. EM in Practice Initial parameters: • Random parameters setting • “Best” guess from other source Stopping criteria: • Small change in likelihood of data • Small change in parameter values Avoiding bad local maxima: • Multiple restarts • Early “pruning” of unpromising ones

16. Sequence Comparison using HMM We now use HMM to extend such log-odds scoring functions to alignments which may contain gaps (indels). Recall: We used log-odds scoring functions for gapless alignments, as follows: s(a,b)=log(pab / qaqb ), where pab and qaare the probabilities of the “Match” and “Random” models.

17. … … (C,-) (G,T) M X Sequence alignment using HMM • Each “output symbol” of the HMM is an aligned pair of two letters, or of a letter and a gap. • Example: Insertion of a first gap in this model: We still need to assign transition/emission probabilities 17

18. Need to define the hidden states, and the transition and emission probabilities, which define the probability of each aligned pair of sequences. • Given two input sequences, we look for an alignment of these sequences of maximum probability.

19. Hidden states and emitted symbols “Hidden” States • Match. • Insertion in x. • insertion in y. Symbols emitted • Match: {(a,b)| a,b in∑ }. • Insertion in x: {(a,-)| a in ∑ }. • Insertion in y: {(-,a)| a in ∑ }.

20. M X Y 1-2δ δ δ M 1- ε ε 0 X 1- ε 0 ε Y Transitions and Emission Probabilities Emission Probabilities • Match: (a,b) with pab – only from M states • Insertion in x: (a,-) with qa – only from X state • Insertion in y: (-,a).with qa - only from Y state. (Note that the hidden states can be reconstructed from the alignment.) • Transitions probabilities • (note the forbidden ones). • δ= probability for 1st gap • ε = probability for tailing gap.

21. Scoring alignments For each pair of sequences x (of length m)and y (of length n), there are many alignments of x and y, each corresponds to a different state-path (the lengths of the paths are between max{m,n} and m+n). Given the transmission and emission probabilities, each alignment has a defined score – the product of the corresponding probabilities. An alignment is “most probable” if it maximizes this score.

22. Finding the most probable alignment Let vM(i,j) be the probability of the most probable alignment of x(1..i) and y(1..j), which ends with a match. Similarly, vX(i,j) and vY(i,j), the probabilities of the most probable alignment of x(1..i) and y(1..j), which ends with an insertion to x or y. Then using a recursive argument, we get:

23. Most probable alignment By similar argument for vX(i,j) and vY(i,j), the probabilities of the most probable alignment of x(1..i) and y(1..j), which ends with an insertion to x or y, are:

24. The Probability Space Different alignments of x and y may have different lengths, so the probability space which we used earlier, for HMM of a fixed length L, is not applicable to this alignment HMM model. However, there is a probability space which contains all infinite sequence alignments (finite alignments are compound events in this model). The algorithm of the previous slides compute the correct probability of each alignment in this probability space. Another approach is to define a probability space which contains all alignments of finite length. In the following we adapt our algorithm to this model.

25. Adding termination probabilities Probability space for all finite alignments is obtained by adding an END state, which denotes the end of the alignment Each other state has a transition probability τ to the END state. This results in an expected sequence length of 1/ τ. The last transition in each alignment is to the END state, with probability τ

26. The log-odds scoring function • We wish to compare the “model” alignment score to the “random” alignment score. • For gapless alignments we used the log-odds ratio: s(a,b) = log (pab / qaqb). • To adapt this for the HMM model, we need to model random sequence by HMM, with end state.

27. scoring function for random model The transition probabilities for the random model, with termination probability η: (x is the start state) The emission probability for a is qa. Thus the probability of x (of length n) and y (of length m) being random is: And the corresponding score is:

28. Markov Matrices for “Random” and “Model” “Random” “Model”

29. Combining models in the log-odds scoring function In order to compare the M score to the R score of sequences x and y, we can find an optimal M score, and then subtract from it the R score. This is insufficient when we look for local alignments, where the optimal substrings in the alignment are not known in advance. A better way: • Define a log-odds scoring function which keeps track of the difference Match-Random scores of the partial strings during the alignment. • At the end add to the score (logτ – 2logη) to compensate for the end transitions in both models. We get the following:

30. The log-odds scoring function (assuming that letters at insertions/deletions are selected by the random model) And at the end add to the score (logτ – 2logη).

31. The log-odds scoring function Another way, with uniform scoring for the M state (Durbin et. al., Chapter 4.1): Define scoring function s with penalties d and e for a first gap and a tailing gap, resp. Then modify the algorithm to correct for extra prepayment, as follows:

32. Log-odds alignment algorithm Initialization: VM(0,0)=logτ - 2logη. Termination: V = max{ VM(m,n), VX(m,n)+c, VY(m,n)+c} Where c= log (1-2δ-τ) – log(1-ε-τ)

33. Total probability of x and y Rather then computing the probability of the most probable alignment, we look for the total probability that x and y are related by our model. Let fM(i,j) be the sum of the probabilities of all alignments of x(1..i) and y(1..j), which end with a match. Similarly, fX(i,j) and fY(i,j) are the sum of these alignments which end with insertion to x (y resp.). A “forward” type algorithm for computing these functions. Initialization: fM(0,0)=1, fX(0,0)= fY(0,0)=0 (we start from M, but we could select other initial state).

34. Total probability of x and y (cont.) The total probability of all alignments is: P(x,y|model)= fM[m,n] + fX[m,n] + fY[m,n]