17

Hidden Markov Models

instructor: Ross A. Lippert

http://www-math.mit.edu/~lippert/18.417/

Announcements:

Start chapter 11
Problem set 6 available

The CG island phenomenon

Nucleotide frequencies in the human genome

a	c	g	t
29.5	20.4	20.5	29.6

aa	ac	ag	at	ca	cc	cg	ct	ga	gc	gg	gt	ta	tc	tg	tt
0.0978	0.0503	0.0699	0.0773	0.0725	0.0521	0.0098	0.07	0.0593	0.0426	0.0521	0.0505	0.0657	0.0594	0.0727	0.098
0.0872	0.0604	0.0604	0.0873	0.0604	0.0418	0.0418	0.0605	0.0604	0.0418	0.0418	0.0605	0.0873	0.0605	0.0605	0.0875

Explanation: CG frequently mutates into TG

requires methylation of the CG
this is suppressed around genes

How can we definitively locate high CG regions?

Numerous heuristics

look for 400bps with >50% CG
look for 50bps with >75% CG

None of these satisfying

Investigate two models of sequence

Normal: CG occurs infrequently and everything else random
CG rich: CG occurs more often

What can help us?

Previously identified CG islands
Known junk sequence

Simplification: two coins Bernoulli sequence

Suppose we have two coins, differently biased

coin used	X	X	X	Y	Y	Y	X	X
result	1	1	0	0	0	0	1	0

Probability problem 1: likelihood of result given hidden information

input: sequence of hidden states, s1,...sn and observations x1,...,xn
output: probability of the x's

What if we don't have the s's, but just have the x's?

A probabilistic model for the s's

Bernoulli just reduces to more Bernoulli
Markov Chain incorporates limited dependency

Markov Models

Markov models generate sequences over an alphabet

Alternative description/notation

Transition matrix T =

P(F->F)	P(F->B)
P(B->F)	P(B->B)

Emission matrix E =

P(F~>H)	P(F~>T)
P(B~>H)	P(B~>T)

Starting vector pi =

P(0->F)

P(0->B)

Using Hidden Markov Models for sequence scoring

Probability problem 2:

input: T and E for HMM and observations x1,...,xn
output: probability of the x's

Formally reduces to a sum of problem 1 over all possible s strings, weighted by the Markov model

This sum is very big, but can be cast in terms of matrix multiplies

pi D(x1) T D(x2) T ... D(xn) 1T = transition matrix, D(x) = diag(E(:,x)), 1 is a vector of 1s

Obtaining conditional probabilities with HMMs

Two traditionally important auxillary quantities

f(k,s) = P(x1...xk,sk=s) = pi D(x1) T ... T D(xk)
b(k,s) = P(x{k+1}...xn,sk=s) = T D(x{k+1}) ... T D(xn) 1

Can now calculate P(sk=s|x1...xn)

Can we use P(sk=s|x1...xn) to infer a sequence of states?

does not produce a consistent sequence

Inferring best state sequence

Maximum likelihood problem 1:

input: HMM with T and E with observations x1...xn
output: s1...sn such that P(s1...sn|x1...xn)

It is really a tropicalized version of P(x1...xn)

(+,*) --> (min,+), prob --> -log(prob)
sum over path products --> min over path sums
ft(k,s) = pi*D(x1)*T*...*T*D(xk)
bt(k,s) = T*D(x{k+1})*...*T*D(xn)*1
sk = argmin f(k,s)+b(k,s)

This is called `the Viterbi algorithm'

Training an HMM

We can score sequences, we can tag states, but where do HMMs come from?

Maximum likelihood problem 2:

input: a (set of) sequence x1...xn and a parameter #s
output: the #s state HMM with max P(x1...xn)

Nonlinear optimization problem with constraints: M >= 0, E >= 0, M 1 = 1, E 1 = 1

P(x1...xn) is a polynomial in M,E
polynomial optimization is NP complete

Expectation maximization

Very commonly used in machine learning

An outgrowth of general considerations of gradient search for constrained polynomials

a^ = avg # of times a is used

P(s~>x) <-- (1/C) SUM_{xk=x} f(k,s) b(k,s)
P(s->s') <-- (1/C) SUM_{k} f(k,s) P(s->s') P(s'->xk) b(k+1,s')

Baum-Welch EM algorithm

Converges to a local optimum

GENSCAN

leading HMM for gene discovery

Global alignment

The maximum likelihood path in this HMM is equivalent to global alignment with affine gap penalties

I = insertion states
D = deletion states
M = matching/mismatching states

Can be generalized to sequence profiles