5481: CSCI Computational Techniques in Genomics (Fall 2008)
Hidden Markov Models for Biological Sequence Analysis
The slides are modified from lecture notes by Serafim Batzoglou and bioalgorithms.info
When the right answer is unknown
We dont know the true Akl, Ek(b) Idea:
We estimate our best guess on what Akl, Ek(b) are We update the parameters of the model, based on our guess We repeat
The Baum-Welch Algorithm
Initialization: Pick the best-guess for model parameters (or arbitrary) Iteration:
1. 2. 3. 4. 5.
Forward Backward Calculate Calculate new model parameters Calculate new log-likelihood
Akl, Ek(b), given CURRENT NEW : akl, ek(b) P(x | NEW)
GUARANTEED TO BE HIGHER BY EXPECTATION-MAXIMIZATION Until P(x | ) does not change much
Alternative: Viterbi Training
Initialization: Iteration:
1. 2. 3.
Same
Perform Viterbi, to find * Calculate Akl, Ek(b) according to * + pseudocounts Calculate the new parameters akl, ek(b)
Until convergence Notes:
Not guaranteed to increase P(x | ) Guaranteed to increase P(x, | , *) In general, worse performance than Baum-Welch
Higher-order HMMs
How do we model memory larger than one time point? P(i+1 = l | i = k) akl P(i+1 = l | i = k, i -1 = j) ajkl A second order HMM with K states is equivalent to a first order HMM with K2 states
aHHT aHT(prev = H) aHT(prev = T)
state HH
aHTH aTHH
state HT
aTHT aTTH aHTT
state H
state T
aTH(prev = H) aTH(prev = T)
state TH
state TT
Solution 1: Chain several states
p 1-p X Y q
X
X
1-q Disadvantage: Still very inflexible lX = C + geometric with mean 1/(1-p)
Solution 2: Negative binomial distribution
p 1p X(1)
Duration in X: m turns, where
p 1p X(2)
p 1p X(n) Y
During first m 1 turns, exactly n 1 arrows to next state are followed During mth turn, an arrow to next state is followed
m1 P(lX = m) = n1 (1 p)n-1+1p(m-1)-(n-1) =
m1 n 1 (1 p)npm-n
HMM interpretation of an alignment
An alignment is a hypothesis that the two sequences are related by evolution Goal: Produce the most likely alignment Assert the likelihood that the sequences are indeed related
A Pair HMM for alignments
1 2 M P(xi, yj) I P(xi) optional J P(yj)
This model generates two sequences simultaneously Match/Mismatch state M: P(x, y) reflects substitution frequencies between pairs of amino acids
1
1
Insertion states I, J: P(x), P(y) reflect frequencies of each amino acid : set so that 1/2 is avg. length before next match : set so that 1/(1 ) is avg. length of a gap
A Pair HMM for unaligned sequences
1 I P(xi) 1 J P(yj)
Two sequences are independently generated from one another P(x, y | R) = P(x1)P(xm) P(y1)P(yn) =
i P(xi) j P(yj)
To compare ALIGNMENT vs. RANDOM hypothesis 1 2
1 I P(xi) M P(xi, yj) J P(yj) 1
1 I
P(xi)
1 J
P(yj)
To compare ALIGNMENT vs. RANDOM hypothesis
Idea: We will divide alignment score by the random score, and take logarithms Let P(xi, yj) s(xi, yj) = log + log (1 2) P(xi) P(yj) (1 ) P(xi) = log (1 2) P(xi) P(xi) = log P(xi)
=Defn substitution score
d
=Defn gap initiation penalty
e
=Defn gap extension penalty
The meaning of alignment scores
The Viterbi algorithm for Pair HMMs corresponds exactly to global alignment DP with affine gaps
VM(i, j) = VI(i, j) = VJ(i, j) = max { VM(i 1, j 1), VI( i 1, j 1) d, Vj( i 1, j 1) } + s(xi, yj) max { VM(i 1, j) d, VI( i 1, j) e } max { VM(i , j - 1) d, Vj( I , j - 1) e }
s(.,.)
~how often a pair of letters substitute one another 1/mean length of next gap 1/mean arrival time of next gap
CSCI 5481: Computational Techniques in Genomics (Fall 2008)
Multiple Sequence Alignment
Rui Kuang
Department of Computer Science and Engineering University of Minnesota (The slides are modified from lecture notes by bioalgorithms.info)
Multiple Alignment versus Pairwise Alignment
Up until now we have only tried to align two sequences. What about more than two? And what for? A faint similarity between two sequences becomes signicant if present in many Multiple alignments can reveal subtle similarities that pairwise alignments do not reveal
Generalizing the Notion of Pairwise Alignment
Alignment of 2 sequences is represented as a 2-row matrix In a similar way, we represent alignment of 3 sequences as a 3-row matrix
A T _ G C G _ A _ C G T _ A A T C A C _ A
Score: more conserved columns, better = alignment
Alignments Paths in
A -T G C
Align 3 sequences: ATGC, AATC,ATGC
A
A
T
--
C
--
A
T
G
C
Alignment Paths
0 1 A 1 -2 T 3 G 4 C
x coordinate
A
A
T
--
C
--
A
T
G
C
Alignment Paths
0 1 A 0 1 A 1 -2 A 2 T 3 T 3 G 3 -4 C 4 C
Align the following 3 sequences: ATGC, AATC,ATGC
x coordinate y coordinate
--
A
T
G
C
Alignment Paths
0 1 A 0 1 A 0 0 -1 -2 A 1 A 2 T 3 T 2 T 3 G 3 -3 G 4 C 4 C 4 C
x coordinate y coordinate z coordinate
Resulting path in (x,y,z) space:
(0,0,0)(1,1,0)(1,2,1) (2,3,2) (3,3,3) (4,4,4)
Aligning Three Sequences
Same strategy as aligning two sequences Use a 3-D Manhattan Cube, with each axis representing a sequence to align For global alignments, go from source to sink
source
sink
2-D vs 3-D Alignment Grid
V W
2-D edit graph
3-D edit graph
2-D cell versus 2-D Alignment Cell
In 2-D, 3 edges in each unit square In 3-D, 7 edges in each unit cube
Architecture of 3-D Alignment Cell
(i-1,j-1,k1) (i-1,j,k-1) (i-1,j-1,k) (i-1,j,k)
(i,j-1,k1)
(i,j,k-1) (i,j-1,k) (i,j,k)
Multiple Alignment: Dynamic Programming
si,j,k = max
si-1,j-1,k-1 + (vi, wj, uk) si-1,j-1,k + (vi, wj, _ ) si-1,j,k-1 + (vi, _, uk) si,j-1,k-1 + (_, wj, uk) si-1,j,k + (vi, _ , _) si,j-1,k + (_, wj, _) si,j,k-1 + (_, _, uk)
cube diagonal: no indels face diagonal: one indel edge diagonal: two indels
(x, y, z) is an entry in the 3-D scoring matrix
Multiple Alignment: Running Time
For
3 sequences of length n, the run time is 7n3; O(n3) k sequences, build a k-dimensional Manhattan, with run time (2k-1)(nk); O (2knk) dynamic programming approach for alignment between two sequences is easily extended to k sequences but it is impractical due to exponential running time
For
Conclusion:
Multiple Alignment Induces Pairwise Alignments
Every multiple alignment induces pairwise alignments
x: y: z: Induces: x: ACGCGG-C; y: ACGC-GAC; x: AC-GCGG-C; z: GCCGC-GAG; y: AC-GCGAG z: GCCGCGAG AC-GCGG-C AC-GC-GAG GCCGC-GAG
Reverse Problem: Constructing Multiple Alignment from Pairwise Alignments
Given 3 arbitrary pairwise alignments: x: ACGCTGG-C; y: ACGC--GAC; x: AC-GCTGG-C; z: GCCGCA-GAG; y: -ACGC-GAG z: GCCGCAGAG
can we construct a multiple alignment that induces them?
Reverse Problem: Constructing Multiple Alignment from Pairwise Alignments
Given 3 arbitrary pairwise alignments: x: ACGCTGG-C; y: ACGC--GAC; x: AC-GCTGG-C; z: GCCGCA-GAG; y: -ACGC-GAG z: GCCGCAGAG
can we construct a multiple alignment that induces them? y: AC-GC-GAG NOT ALWAYS
z: GCCGCAGAG
Pairwise alignments may be inconsistent
Combining Optimal Pairwise Alignments...