Markov Chains Lecture #5

1. Markov ChainsLecture #5 Background Readings: Durbin et. al. Section 3.1 Prepared by Shlomo Moran, based on Danny Geiger’s and Nir Friedman’s .

2. Dependencies along the genome • In previous classes we assumed every letter in a sequence is sampled randomly from some distribution q() over the alpha bet {A,C,T,G}. • This is too restrictive for true genomes. • There are special subsequences in the genome, like TATA within the regulatory area, upstream a gene. • The pattern CG is less common than expected for random sampling. • We model such dependencies by Markov chains and hidden Markov model, which we define next.

3. Finite Markov Chain • An integer time stochastic process, consisting of a domainD of m>1 states {s1,…,sm} and • An m dimensional initial distribution vector ( p(s1),.., p(sm)). • Anm×mtransition probabilities matrix M= (msisj) • For example, D can be the letters {A, C, T, G}, p(A) the probability of A to be the 1st letter in a sequence, and mAG the probability that G follows A in a sequence.

4. X1 X2 Xn-1 Xn • For each integer n, a Markov Chain assigns probability to sequences (x1…xn) inDn (i.e, xiD) as follows: Markov Chain (cont.)

5. X1 X2 Xn-1 Xn Markov Chain (cont.) Similarly, each Xi is a probability distributions over D, which is determined by the initial distribution (p1,..,pn) and the transition matrix M. There is a rich theory which studies the properties of such “Markov sequences”(X1,…, Xi ,…). A bit of this theory is presented next.

6. A B C D 0.95 0 0.05 0 A 0.2 0.5 0 0.3 B 0 0.2 0 0.8 C 0 0 1 0 D Matrix Representation The transition probabilities Matrix M =(mst) M is a stochastic Matrix: The initial distribution vector (u1…um)defines thedistributionof X1 (p(X1=si)=ui). After one move, the distribution is changed to X2 = X1M

7. A B C D 0.95 0 0.05 0 A 0.2 0.5 0 0.3 B 0 0.2 0 0.8 C 0 0 1 0 D Matrix Representation Example: if X1=(0, 1, 0, 0) then X2=(0.2, 0.5, 0, 0.3) And if X1=(0, 0, 0.5, 0.5) then X2=(0, 0.1, 0.5, 0.4). The i-th distribution is Xi= X1Mi-1

8. 0.95 A B C D 0.95 0 0.05 0 0.2 0.5 A A B 0.2 0.5 0 0.3 B 0.2 0.3 0.05 0.8 0 0.2 0 0.8 D C C 1 0 0 1 0 D Representation of a Markov Chain as a Digraph Each directed edge AB is associated with the positive transition probability from A to B.

9. B A D C Properties of Markov Chain states States of Markov chains are classified by the digraph representation (omitting the actual probability values) A, C and D are recurrent states: they are in strongly connected components which are sinks in the graph. B is not recurrent – it is a transient state Alternative definitions: A state s is recurrent if it can be reached from any state reachable from s; otherwise it is transient.

10. A B D C Another example of Recurrent and Transient States A and B are transient states, C and D are recurrent states. Once the process moves from B to D, it will never come back.

11. E A B D C Irreducible Markov Chains A Markov Chain is irreducible if the corresponding graph is strongly connected (and thus all its states are recurrent). A B D C

12. E A B D C Periodic States A state s has a period k if k is the GCD of the lengths of all the cycles that pass via s. (in the shown graph the period of A is 2). Exercise: All the states in the same strongly connected component have the same period A Markov Chain is periodic if all the states in it have a period k >1. It is aperiodic otherwise.

13. Ergodic Markov Chains A B • A Markov chain is ergodic if : • the corresponding graph is strongly connected. • It is not peridoic D C Ergodic Markov Chains are important since they guarantee the corresponding Markovian process converges to a unique distribution, in which each state has a strictly positive probability.

14. Stationary Distributions for Markov Chains Let M be a Markov Chain of m states, and let V = (v1,…,vm) be a probability distribution over the m states V = (v1,…,vm) is stationary distribution for M if VM=V. (ie, if one step of the process does not change the distribution). V is a stationary distribution V is a non-negative left (row) Eigenvector of M with Eigenvalue 1.

15. Stationary Distributions for a Markov Chain M Exercise: A stochastic matrix always has a real left Eigenvector with Eigenvalue 1 (hint: show that a stochastic matrix has a right Eigenvector with Eigenvalue 1. Note that the left Eigenvalues of a Matrix are the same as the right Eigenvlues). [It can be shown that the above Eigenvector V can be non-negative. Hence each Markov Chain has a stationary distribution.]

16. “Good” Markov chains • A Markov Chains is good if the distributions Xi , as i∞: • converge to a unique distribution, independent of the initial distribution. • (2) In that unique distribution, each state has a positive probability. • The Fundamental Theorem of Finite Markov Chains:A Markov Chain is good  the corresponding graph is ergodic. • We will prove the  part, by showing that non-ergodic Markov Chains are not good.

17. Examples of “Bad” Markov Chains • A Markov chains is “bad” if either : • It does not converge to a unique distribution. • It does converge to u.d., but some states in this distribution have zero probability.

18. B A D C Bad case 1: Mutual Unreachabaility • Consider two initial distributions: • p(X1=A)=1 (p(X1 = x)=0 if x≠A). • p(X1= C) = 1 In case a), the sequence will stay atA forever. In case b), it will stay in {C,D} for ever. Fact 1: If G has two states which are unreachable from each other, then {Xi} cannot converge to a distribution which is independent on the initial distribution.

19. A B D C Bad case 2: Transient States Once the process moves from B to D, it will never come back.

20. A B D C Bad case 2: Transient States Fact 2: For each initial distribution, with probability 1 a transient state will be visited only a finite number of times. X Proof: Let A be a transient state, and let X be the set of states from which A is unreachable. It is enough to show that, starting from any state, with probability 1 a state in X is reached after a finite number of steps (Exercise: complete the proof)

21. Corollary: A good Markov Chain is irreducible

22. E A B D C Bad case 3: Periodic Markov Chains Recall: A Markov Chain is periodic if all the states in it have a period k >1. The above chain has period 2. In the above chain, consider the initial distribution p(B)=1. Then states {B, C} are visited (with positive probability) only in odd steps, and states {A, D, E} are visited in only even steps.

23. E A B D C Bad case 3: Periodic States Fact 3: In a periodic Markov Chain (of period k >1) there are initial distributions under which the states are visited in a periodic manner. Under such initial distributions Xidoes not converge as i∞. Corollary: A good Markov Chain is not periodic

24. The Fundamental Theorem of Finite Markov Chains: We have proved that non-ergodic Markov Chains are not good. A proof of the opposite direction below is based on Perron-Frobenius theory on nonnegative matrices, and is beyond the scope of this course: • If a Markov Chain is ergodic, then • It has a unique stationary distribution vector V > 0, which is an Eigenvector of the transition matrix. • For any initial distribution, the distributions Xi , as i∞, converges to V.

25. Use of Markov Chains in Genome search: Modeling CpG Islands In human genomes the pair CG often transforms to (methyl-C) G which often transforms to TG. Hence the pair CG appears less than expected from what is expected from the independent frequencies of C and G alone. Due to biological reasons, this process is sometimes suppressed in short stretches of genomes such as in the start regions of many genes. These areas are called CpG islands (p denotes “pair”).

26. Example: CpG Island (Cont.) We consider two questions: Question 1: Given a short stretch of genomic data, does it come from a CpG island ? Question 2: Given a long piece of genomic data, does it contain CpG islands in it, where, what length ? We “solve” the first question by modeling strings with and without CpG islands as Markov Chains over the same states {A,C,G,T} but different transition probabilities:

27. Example: CpG Island (Cont.) The “+” model: Use transition matrix M+ = (m+st), Where: m+st = (the probability that t follows s in a CpG island) The “-” model: Use transition matrix M- = (m-st), Where: m-st = (the probability that t follows s in a non CpG island)

28. Example: CpG Island (Cont.) With this model, to solve Question 1 we need to decide whether a given short sequence of letters is more likely to come from the “+” model or from the “–” model. This is done by using the definitions of Markov Chain, in which the parameters are determined by known data and the log odds-ratio test.

29. Question 1: Using two Markov chains M+ (For CpG islands): We need to specify p+(xi | xi-1) where + stands for CpG Island. From Durbin et al we have: Xi Xi-1 (Recall: rows must add up to one; columns need not.)

30. Question 1: Using two Markov chains M- (For non-CpG islands): …and for p-(xi | xi-1) (where “-” stands for Non CpG island) we have: Xi Xi-1

31. X1 X2 XL-1 XL Discriminating between the two models Given a string x=(x1….xL), now compute the ratio If RATIO>1, CpG island is more likely. Actually – the log of this ratio is computed: Note: p+(x1|x0) is defined for convenience as p+(x1). p-(x1|x0) is defined for convenience as p-(x1).

32. Log Odds-Ratio test Taking logarithm yields If logQ > 0, then + is more likely (CpG island). If logQ < 0, then - is more likely (non-CpG island).

33. Where do the parameters (transition- probabilities) come from ? Learning from complete data, namely, when the label is given and every xi is measured: Source: A collection of sequences from CpG islands, and a collection of sequences from non-CpG islands. Input: Tuples of the form (x1, …, xL, h), where h is + or - Output: Maximum Likelihood parameters (MLE) Count all pairs (Xi=a, Xi-1=b) with label +, and with label -, say the numbers are Nba,+ and Nba,- .

34. X1 X2 XL-1 XL Where Nbx,+ is the number of times letter x appears after letter b in CpG islands in the dataset. Maximum Likelihood Estimate (MLE) of the parameters (using labeled data) The needed parameters are: P+(x1), p+ (xi | xi-1), p-(x1), p-(xi | xi-1) The ML estimates are given by: Where Nx,+ is the number of times letter x appear in CpG islands in the dataset. Using MLE is justified when we have a large sample. The numbers appearing in our tables (taken from Durbin et al p. 50) are based on 60,000 nucleotides.

35. Question 2: Finding CpG Islands Given a long genomic string with possible CpG Islands, we define a Markov Chain over 8 states, all interconnected (hence it is ergodic): The problem is that we don’t know the sequence of states which are traversed, but just the sequence of letters. A+ C+ G+ T+ A- C- G- T- Therefore we use here Hidden Markov Model

36. M M M M S1 S2 SL-1 SL T T T T x1 x2 XL-1 xL Hidden Markov Model A Markov chain (s1,…,sL): and for each state s and a symbol x we have p(Xi=x|Si=s) Application in communication: message sent is (s1,…,sm) but we receive (x1,…,xm) . Compute what is the most likely message sent. Application in speech recognition: word said is (s1,…,sm) but we recorded (x1,…,xm) . Compute what is the most likely word said.

37. M M M M S1 S2 SL-1 SL T T T T x1 x2 XL-1 xL Hidden Markov Model Notations: Markov Chain transition probabilities: p(Si+1= t|Si= s) = ast Emission probabilities: p(Xi = b| Si = s) = es(b) For Markov Chains we know: What is p(s,x) = p(s1,…,sL;x1,…,xL)?

38. M M M M S1 S2 SL-1 SL T T T T x1 x2 XL-1 xL Hidden Markov Model p(Xi = b| Si = s) = es(b), means that the probability of xidepends only on the probability of si. Formally, this is equivalent to the conditional independence assumption: p(Xi=xi|x1,..,xi-1,xi+1,..,xL,s1,..,si,..,sL) = esi(xi) Thus