Pair HMMs and Profile HMMs COMPSCI 260 Spring 2016 HMM - - PowerPoint PPT Presentation

Pair ¡HMMs ¡and ¡Profile ¡HMMs

¡

COMPSCI ¡260 ¡– ¡Spring ¡2016 ¡

HMM ¡

q 0

100% 80% 15% 30% 70% 5% R=0% Y = 100%

q1

Y=0% R = 100%

q2 M1=({q0,q1,q2},{Y,R},Pt,Pe) Pt={(q0,q1,1), (q1,q1,0.8), (q1,q2,0.15), (q1,q0,0.05), (q2,q2,0.7), (q2,q1,0.3)} Pe={(q1,Y,1), (q1,R,0), (q2,Y,0), (q2,R,1)}

An Example

Three ¡views ¡of ¡an ¡HMM ¡ ¡

QuesAons ¡we ¡can ¡address ¡with ¡an ¡HMM ¡

INPUT ¡ The ¡HMM ¡model ¡M: ¡Q={q0, q1, ... , qm}; Pt (qj | qi); Pe (sj | qi) ¡ A ¡sequence ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡and ¡a ¡path ¡ ¡ OUTPUT ¡ What ¡is ¡the ¡probability ¡of ¡generaAng ¡sequence ¡S ¡from ¡path ¡ϕ according ¡to ¡the ¡model ¡M? ¡P(S|ϕ,M)

q 0

100% 80% 15% 30% 70% 5%

A=25% C=25% G=25% T=25% A=10% C=40% G=10% T=40%

CGATATTCGATTCTACGCGCGTATACTAGCTTATCTGATC 011111112222222111111222211111112222111110

q1 q2

P(S |φ) = P

e(xi | yi+1) i=0 L−1

∏

emission prob.

φ = y0y1...yL+1 S = x0x1...xL−1

OUTPUT ¡ What ¡is ¡the ¡joint ¡probability ¡of ¡sequence ¡S ¡and ¡path ¡ϕ ¡according ¡ to ¡the ¡model ¡M? ¡P(S,ϕ|M)

QuesAons ¡we ¡can ¡address ¡with ¡an ¡HMM ¡

P(S |φ) = P

e(xi | yi+1) i=0 L−1

∏

emission prob.

P(φ) = P

t(yi+1 | yi) i=0 L

∏

transition prob.

P(S,φ) = P(S φ)P(φ)

CGATATTCGATTCTACGCGCGTATACTAGCTTATCTGATC 011111112222222111111222211111112222111110 INPUT ¡ The ¡HMM ¡model ¡M: ¡Q={q0, q1, ... , qm}; Pt (qj | qi); Pe (sj | qi) ¡ A ¡sequence ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡and ¡a ¡path ¡ ¡

φ = y0y1...yL+1 S = x0x1...xL−1

QuesAons ¡we ¡can ¡address ¡with ¡an ¡HMM ¡

INPUT ¡ The ¡HMM ¡model ¡M: ¡Q={q0, q1, ... , qm}; Pt (qj | qi); Pe (sj | qi) ¡ A ¡sequence ¡S: ¡ OUTPUT ¡ What ¡is ¡the ¡most ¡probable ¡path ¡for ¡generaAng ¡sequence ¡S ¡ according ¡to ¡the ¡model ¡M ? -­‑ ¡ ¡ ¡“Decoding” ¡

q 0

100% 80% 15% 30% 70% 5%

A=25% C=25% G=25% T=25% A=10% C=40% G=10% T=40%

CGATATTCGATTCTACGCGCGTATACTAGCTTATCTGATC

q1 q2

φmax = argmax

φ

P(φ S, M)

φmax = argmax

φ

P(φ S) = argmax

φ

P(φ,S) P(S) = argmax

φ

P(φ,S) = argmax

φ

P(S φ)P(φ)

P(φ) = P

t(yi+1 | yi) i=0 L

∏

P(S |φ) = P

e(xi | yi+1) i=0 L−1

∏

φmax = argmax

φ

P

t(q0 yL)

P

e(xi yi+1)P t(yi+1 yi) i=0 L−1

∏

emission prob. transition prob.

φ = y0y1...yL+1 S = x0x1...xL−1

“Decoding” ¡with ¡an ¡HMM ¡– ¡Viterbi ¡decoding ¡

• Viterbi ¡gives ¡us ¡two ¡things: ¡

– The ¡“best” ¡parse: ¡ϕmax = argmaxϕ P(ϕ|S) – The ¡joint ¡probability: ¡P(ϕmax, S)

“Decoding” ¡with ¡an ¡HMM ¡

• What ¡if ¡we ¡are ¡interested ¡in ¡the ¡state ¡that ¡generated ¡a ¡

parAcular ¡character? ¡

P(S)

• What ¡if ¡we ¡are ¡interested ¡in ¡the ¡marginal ¡probability ¡of ¡

emiRng ¡S, ¡regardless ¡of ¡the ¡path? ¡

P(yk = qi S) = P(S,yk = qi) P(S)

• We ¡can ¡compute ¡these ¡using ¡the ¡Forward ¡and ¡Backward ¡

algorithms, ¡and ¡“posterior” ¡decoding ¡

“Decoding” ¡with ¡an ¡HMM ¡-­‑ ¡Posterior ¡decoding ¡

P(yk = qi S) = P(S,yk = qi) P(S) = F(i,k)B(i,k) P(S)

“Posterior” decoding:

Training ¡an ¡HMM ¡with ¡labeled ¡sequences: ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡given ¡ ¡

CGATATTCGATTCTACGCGCGTATACTAGCTTATCTGATC 011111112222222111111222211111112222111110

to state 1 2 from state 0 0 (0%) 1 (100%)

0 (0%)

1 1 (4%) 21 (84%) 3 (12%) 2 0 (0%) 3 (20%) 12 (80%) symbol A C G T in state

1

6 (24%) 7 (28%) 5 (20%) 7 (28%) 2 3 (20%) 3 (20%) 2 (13%) 7 (47%)

∑

− =

=

1 | | , , , Q h h i j i j i

A A a

∑

− Σ =

=

1 | | , , , h h i k i k i

E E e

transitions emissions

S,φ

OUTPUT: ¡ ¡The ¡parameters ¡of ¡the ¡HMM: Pt (qj | qi); Pe (sj | qi) INPUT: ¡ ¡A ¡set ¡of ¡sequences ¡S ¡generated ¡by ¡the ¡HMM; ¡Q={q0, q1, ... , qm} ¡ Two ¡Solu5ons: ¡

• 1. ¡Viterbi ¡Training: ¡Start ¡with ¡random ¡HMM ¡parameters. ¡Use ¡Viterbi ¡to ¡

find ¡the ¡most ¡probable ¡path ¡for ¡each ¡training ¡sequence, ¡and ¡then ¡label ¡ the ¡sequence ¡with ¡that ¡path. ¡ ¡Use ¡labeled ¡sequence ¡training ¡on ¡the ¡ resulAng ¡set ¡of ¡sequences ¡and ¡paths. ¡Iterate ¡unAl ¡Viterbi ¡paths ¡do ¡not ¡

• change. ¡
• 2. ¡Baum-­‑Welch ¡Training: ¡Start ¡with ¡random ¡HMM ¡parameters. ¡Use ¡

posterior ¡decoding ¡to ¡compute ¡the ¡‘forward’ ¡and ¡‘backward’ ¡

• probabiliAes. ¡Sum ¡over ¡all ¡possible ¡paths ¡(rather ¡than ¡the ¡single ¡most ¡

probable ¡one) ¡to ¡esAmate ¡expected ¡counts ¡ ¡Ai, j and ¡Ei,k; ¡then ¡use ¡the ¡ same ¡formulas ¡as ¡for ¡labeled ¡sequence ¡training ¡on ¡these ¡expected ¡

• counts. ¡Iterate ¡unAl ¡the ¡change ¡in ¡P(S | M) < ε

Training ¡an ¡HMM ¡with ¡unlabeled ¡sequences: ¡only ¡ ¡ ¡ ¡ ¡given ¡ ¡

S

HMMs ¡for ¡gene ¡predicAon ¡

HMMs ¡and ¡sequence ¡alignments ¡– ¡Pair ¡HMMs ¡

A ¡Pair ¡HMM ¡is ¡an ¡HMM ¡which ¡has ¡two ¡output ¡channels ¡rather ¡than ¡

• ne; ¡each ¡state ¡can ¡emit ¡a ¡symbol ¡into ¡one ¡or ¡the ¡other ¡(or ¡both) ¡
• channels. ¡

More ¡general ¡Pair ¡HMMs ¡can ¡have ¡many ¡more ¡states, ¡but ¡those ¡states ¡can ¡ all ¡be ¡classified ¡as ¡inser7on ¡states, ¡dele7on ¡states, ¡or ¡match/mismatch ¡

• states. ¡

HMMs ¡and ¡sequence ¡alignments ¡– ¡Pair ¡HMMs ¡

q0 q0 I

The ¡most ¡probable ¡state ¡path ¡ through ¡the ¡Pair ¡HMM ¡ determines ¡the ¡op5mal ¡

• alignment. ¡

M D

A ¡Simple ¡Pair ¡HMM ¡ ¡

• ------AACGCAGGAGCCTGCAGGTCTGGGCAGCCAGTTAGCGGGCTGCGGGCCCAGGA bosTau2 0 60 + .

CACTCCCAT--------------GGCCCGG--AGCC------CGAGCCGCGCGCCCACAA canFam2 0 60 + . AGCCCTCGCAGAGCCCTGGGAGAGACAGCCTACAGGACTGGACTTGGGGCAGGGAAACAT bosTau2 60 60 + .

• --CCTGGCAGAGCGCCGGGAGCCGCAGCCTCCAGACCCGAGCGCGCAGGCGGCAGAACG canFam2 60 60 + .

TTCAGAGAAAAGATAGGAGATA bosTau2 120 22 + . CGCGGAG---GGGCGGGCGCCA canFam2 120 22 + .

HMMs ¡and ¡sequence ¡alignments ¡– ¡Pair ¡HMMs ¡

Pair ¡HMMs ¡can ¡be ¡used ¡for ¡simultaneous ¡alignment ¡and ¡annotaAon ¡

q0 q0 I I M D M D

A ¡Pair ¡HMM ¡with ¡ FuncAonal ¡States ¡

• ------AACGCAGGAGCCTGCAGGTCTGGGCAGCCAGTTAGCGGGCTGCGGGCCCAGGA bosTau2 0 60 + .

CACTCCCAT--------------GGCCCGG--AGCC------CGAGCCGCGCGCCCACAA canFam2 0 60 + . AGCCCTCGCAGAGCCCTGGGAGAGACAGCCTACAGGACTGGACTTGGGGCAGGGAAACAT bosTau2 60 60 + .

• --CCTGGCAGAGCGCCGGGAGCCGCAGCCTCCAGACCCGAGCGCGCAGGCGGCAGAACG canFam2 60 60 + .

TTCAGAGAAAAGATAGGAGATA bosTau2 120 22 + . CGCGGAG---GGGCGGGCGCCA canFam2 120 22 + .

Generalization: Profile HMMs

Profile ¡HMMs ¡applicaAon: ¡Pfam ¡protein ¡domains ¡

Profile ¡HMMs ¡applicaAon: ¡Pfam ¡protein ¡domains ¡

PosiAon ¡weight ¡matrices ¡(PWMs) ¡(PSSMs) ¡

G A T C T C A T T T

State ¡transiAon ¡diagram ¡

q0 q0

Graphical ¡model ¡

PWMs ¡are ¡a ¡special ¡case ¡of ¡an ¡HMM: ¡

What are the transition probabilities?

Profile ¡HMMs ¡for ¡protein ¡families ¡

• Consider ¡the ¡PWM ¡for ¡a ¡conserved ¡segment ¡of ¡a ¡protein ¡family ¡
• The ¡profile ¡consists ¡of ¡the ¡frequencies ¡of ¡amino ¡acids ¡at ¡each ¡posiAon ¡ ¡

P(L) ¡= ¡2/3 ¡ P(R) ¡= ¡1/3 ¡ ¡

• However, ¡this ¡type ¡of ¡profile ¡does ¡not ¡allow ¡for ¡gaps ¡(inserAons/

deleAons) ¡

Begin ¡ M1 ¡ M2 ¡ M6 ¡ End ¡ M3 ¡ M4 ¡ M5 ¡ P(R) ¡= ¡2/3 ¡ P(L) ¡= ¡1/3 ¡ P(A) ¡= ¡1/3 ¡ P(V) ¡= ¡1/3 ¡ P(R) ¡= ¡1/3 ¡ ¡ P(Y) ¡= ¡1 ¡ P(V) ¡= ¡2/3 ¡ P(A) ¡= ¡1/3 ¡ P(I) ¡= ¡1 ¡

R I Y V R

Profile ¡HMMs ¡for ¡protein ¡families ¡

Begin ¡ M1 ¡ M2 ¡ ML ¡ End ¡ … I1 I2 I0 IL

• Adding ¡HMM ¡states ¡that ¡capture ¡inserAons: ¡
• For ¡any ¡parAcular ¡insert ¡state, ¡we ¡may ¡have ¡different ¡transiAon ¡

probabiliAes ¡for ¡entering ¡it ¡the ¡first ¡Ame ¡vs. ¡staying ¡in ¡the ¡insert ¡state ¡

• What ¡kind ¡of ¡gap ¡penalty ¡model ¡can ¡Profile ¡HMMs ¡capture? ¡
• Affine ¡gap ¡penalty ¡ ¡
• What ¡is ¡the ¡length ¡distribuAon ¡of ¡the ¡inserted ¡segments ¡(number ¡of ¡

elements ¡emiked ¡from ¡an ¡insert ¡state ¡Ik) ¡? ¡

• Geometric ¡distribuAon ¡

Profile ¡HMMs ¡for ¡protein ¡families ¡

Begin ¡ M1 ¡ M2 ¡ ML ¡ End ¡ … M3 ¡ D1 ¡ D2 ¡ D3 ¡ DL ¡

• Adding ¡HMM ¡states ¡that ¡capture ¡deleAons ¡

Begin ¡ M1 ¡ M2 ¡ ML ¡ End ¡ … M3 ¡

• One ¡could ¡model ¡ ¡

deleAons ¡as: ¡

• However, ¡arbitrarily ¡long ¡gaps ¡introduce ¡lots ¡of ¡transiAons ¡in ¡the ¡model ¡
• Instead, ¡we ¡will ¡introduce ¡delete ¡states ¡that ¡do ¡not ¡emit ¡any ¡symbols ¡ ¡

Profile ¡HMMs ¡for ¡protein ¡families ¡

• The ¡full ¡Profile ¡HMM ¡model: ¡

Begin ¡ M1 ¡ M2 ¡ ML ¡ End ¡ … I1 I2 I0 IL D1 ¡ D2 ¡ DL ¡ … …

• The ¡model ¡incorporates ¡states ¡for ¡inser7ons ¡and ¡dele7ons, ¡as ¡well ¡

as ¡L ¡match/mismatch ¡states ¡ ¡

• Have ¡we ¡finished ¡building ¡the ¡model? ¡
• No, ¡we ¡have ¡just ¡given ¡the ¡overall ¡topology ¡of ¡a ¡profile-­‑HMM, ¡but ¡

we ¡sAll ¡need ¡to ¡decide ¡how ¡many ¡states ¡our ¡HMM ¡has ¡(L), ¡what ¡are ¡ the ¡transi7on ¡probabili7es, ¡what ¡are ¡the ¡emission ¡probabili7es. ¡

Profile ¡HMMs ¡for ¡protein ¡families ¡

Begin ¡ M1 ¡ M2 ¡ ML ¡ End ¡ … I1 I2 I0 IL D1 ¡ D2 ¡ DL ¡ … …

• A ¡few ¡important ¡points ¡to ¡keep ¡in ¡mind ¡ ¡

– The ¡delete ¡states ¡emit ¡no ¡characters ¡ ¡ – Emission ¡probabili7es ¡exist ¡only ¡for ¡insert ¡and ¡match/mismatch ¡ states ¡ – Transi7on ¡probabili7es ¡between ¡states ¡determine ¡how ¡likely ¡it ¡ is ¡to ¡insert ¡or ¡delete ¡an ¡element ¡

Profile ¡HMMs ¡for ¡protein ¡families ¡

V G A

• H

A G E Y V

• N

V D E V V E A

• D

V A G H V K G

• D

V Y S

• T

Y E T S F N A

• N

I P K H I A G A D N G A G Y

1 ¡ 2 ¡ 3 ¡ 4 ¡ 5 ¡ 6 ¡ 7 ¡ 8 ¡ Begin ¡ M1 ¡ M2 ¡ M8 ¡ End ¡ … I1 I2 I0 I8

D1 ¡ D2 ¡ D8 ¡

… … M3 ¡ M4 ¡ I3 I4

D3 ¡ D4 ¡

Pfam ¡protein ¡domains ¡

Pfam ¡protein ¡domains ¡

Column height is relative entropy (or information content) at that position: Information content (bits) Letter heights are proportional to frequency Insertion probability

Pfam ¡protein ¡domains ¡

Profile ¡HMMs ¡-­‑ ¡training ¡

Begin ¡ M1 ¡ M2 ¡ ML ¡ End ¡ … I1 I2 I0 IL D1 ¡ D2 ¡ DL ¡ … …

• Given ¡the ¡topology ¡of ¡a ¡profile ¡HMM, ¡what ¡algorithm(s) ¡do ¡we ¡use ¡to ¡

train ¡the ¡model? ¡

• It ¡depends ¡on ¡what ¡kind ¡of ¡training ¡data ¡is ¡available: ¡labeled ¡or ¡unlabeled ¡
• Labeled ¡training ¡data: ¡counAng ¡(maximum ¡likelihood ¡approach) ¡
• Unlabeled ¡training ¡data: ¡the ¡Baum-­‑Welch ¡algorithm ¡

Profile ¡HMMs ¡– ¡training ¡using ¡labeled ¡(?) ¡data ¡

VGA--HAGEY V----NVDEV VEA--DVAGH VKG------D VYS--TYETS FNA--NIPKH IAGADNGAGY

• How ¡do ¡we ¡pick ¡the ¡length ¡of ¡the ¡HMM ¡(i.e., ¡how ¡

many ¡match ¡states ¡do ¡we ¡have ¡in ¡the ¡profile)? ¡

• One ¡heurisAc ¡is ¡to ¡only ¡include ¡those ¡columns ¡

that ¡have ¡amino ¡acids ¡in ¡at ¡least ¡half ¡of ¡the ¡ sequences ¡

• For ¡example, ¡in ¡the ¡above ¡alignment, ¡there ¡would ¡be ¡match ¡states ¡

for ¡each ¡column ¡except ¡for ¡the ¡fourth ¡and ¡fioh ¡columns ¡ ¡

• Another ¡reasonable ¡guess ¡is ¡the ¡mean ¡or ¡median ¡length ¡of ¡all ¡the ¡

input ¡sequences ¡

VGA--HAGEY V----NVDEV VEA--DVAGH VKG------D VYS--TYETS FNA--NIPKH IAGADNGAGY

Profile ¡HMMs ¡– ¡training ¡using ¡labeled ¡(?) ¡data ¡

• How ¡do ¡we ¡pick ¡the ¡length ¡of ¡the ¡HMM ¡(i.e., ¡how ¡

many ¡match ¡states ¡do ¡we ¡have ¡in ¡the ¡profile)? ¡

• One ¡heurisAc ¡is ¡to ¡only ¡include ¡those ¡columns ¡

that ¡have ¡amino ¡acids ¡in ¡at ¡least ¡half ¡of ¡the ¡ sequences ¡

• How ¡do ¡we ¡compute ¡emission ¡probabili7es? ¡
• Counts ¡– ¡similarly ¡to ¡PWMs ¡(we ¡may ¡want ¡to ¡add ¡pseudocounts) ¡ ¡

ei,k = Ei,k Ei,h

h=0 |α|−1

∑

• How ¡do ¡we ¡compute ¡transi7on ¡probabili7es? ¡
• Counts ¡

∑

− =

=

1 | | , , , Q h h i j i j i

A A a

Profile ¡HMMs ¡– ¡training ¡using ¡unlabeled ¡data ¡

Begin ¡ M1 ¡ M2 ¡ ML ¡ End ¡ … I1 I2 I0 IL D1 ¡ D2 ¡ DL ¡

VGAHAGEY VNVDEV VEADVAGH VKGD VYSTYETS FNANIPKH IAGNGAGY

• We ¡can ¡use ¡the ¡Baum-­‑Welch ¡algorithm ¡
• Note, ¡however, ¡that ¡the ¡topology ¡of ¡the ¡HMM ¡is ¡

fixed ¡before ¡learning ¡

• In ¡pracAce, ¡Profile ¡HMMs ¡are ¡typically ¡learned ¡from ¡

aligned ¡sequences ¡

Using ¡Profile ¡HMMs ¡to ¡idenAfy ¡and ¡align ¡new ¡family ¡members ¡

• Given ¡a ¡Profile ¡HMM ¡for ¡a ¡protein ¡family, ¡we ¡can ¡use ¡it ¡to ¡query ¡the ¡

GENBANK ¡or ¡some ¡other ¡database ¡to ¡find ¡new ¡family ¡members ¡ ¡

• Given ¡the ¡parameters ¡of ¡the ¡model ¡(Pt (qj | qi) and Pe (sj | qi)) ¡and ¡a ¡

new ¡sequence ¡S, ¡it ¡is ¡possible ¡to ¡determine ¡the ¡most ¡probable ¡ sequence ¡of ¡states ¡for ¡generaAng ¡S ¡… ¡ ¡… ¡using ¡Viterbi ¡decoding ¡

• By ¡doing ¡that, ¡we ¡automaAcally ¡align ¡the ¡new ¡sequence ¡to ¡the ¡
• ther ¡family ¡members ¡
• It ¡is ¡also ¡possible ¡to ¡determine ¡the ¡origin ¡of ¡a ¡parAcular ¡amino ¡acid ¡

(whether ¡it ¡is ¡from ¡some ¡match ¡column ¡j, ¡or ¡an ¡inserAon) ¡…. ¡ ¡… ¡using ¡posterior ¡decoding ¡

• Using ¡either ¡of ¡the ¡two ¡methods, ¡we ¡can ¡select ¡those ¡sequences ¡S ¡

that ¡have ¡large ¡probability ¡of ¡being ¡generated ¡by ¡the ¡given ¡profile ¡ ¡

• But ¡where ¡does ¡the ¡Profile ¡HMM ¡come ¡from? ¡

¡

Chicken ¡and ¡egg ¡problem ¡… ¡

• Given ¡a ¡mul7ple ¡sequence ¡alignment ¡between ¡members ¡of ¡a ¡

protein ¡family, ¡we ¡can ¡build ¡a ¡Profile ¡HMM ¡for ¡the ¡family ¡

• Given ¡a ¡Profile ¡HMM, ¡we ¡can ¡find ¡the ¡mulAple ¡sequence ¡

alignment ¡between ¡members ¡of ¡the ¡family ¡

We ¡can ¡use ¡an ¡iteraAve ¡approach ¡

• We ¡start ¡with ¡a ¡single ¡sequence ¡S ¡and ¡use ¡BLAST ¡to ¡get ¡more ¡

sequences ¡like ¡S ¡

• We ¡use ¡these ¡sequences ¡to ¡construct ¡a ¡profile ¡representaAon ¡
• We ¡use ¡the ¡profile ¡to ¡search ¡the ¡database ¡again ¡for ¡similar ¡sequences

¡

• The ¡new ¡sequences ¡are ¡used ¡to ¡improve ¡the ¡profile ¡representaAon ¡

and ¡the ¡process ¡is ¡repeated… ¡ ¡

• This ¡is ¡the ¡principle ¡behind ¡PSIBLAST ¡-­‑ ¡Posi5on ¡Specific ¡ ¡

Iterated ¡BLAST ¡

HMMER is a software suite for protein sequence similarity searches using probabilistic methods. …. Methods are available for searching either a single protein sequence, multiple protein sequence alignment or profile HMM against a target sequence database, and for searching a protein sequence against Pfam…

More ¡informaAon…. ¡

http://hmmer.org