Transcript Lecture 5 slides

CS 5263 Bioinformatics
Lecture 5: Affine Gap Penalties
Last lecture
• Local Sequence Alignment
• Bounded Dynamic Programming
• Linear Space Sequence Alignment
The Smith-Waterman algorithm
Initialization: F(0, j) = F(i, 0) = 0
0
F(i – 1, j) – d
Iteration: F(i, j) = max
F(i, j – 1) – d
F(i – 1, j – 1) + (xi, yj)
The Smith-Waterman algorithm
Termination:
1. If we want the best local alignment…
FOPT = maxi,j F(i, j)
2. If we want all local alignments scoring > t
For all i, j find F(i, j) > t, and trace back
Bounded Dynamic Programming
yN ………………………… y1
x1 ………………………… xM
• O(kM) time
• O(kN) memory
k
• O(M+N) memory
• 2MN time
k*
M/2
M/2
Linear-space alignment
N-k*
Homework Problem 5 hints
Sequence 2
50
150
200
250
300
50
100
Sequence 1
Dot matrix for visualizing seq similarities
• Seq1: x[1..m]
• Seq2: y[1..n]
100
150
200
A(i, j) = 1 if (xi, yj) = 1
250
Sequence 2
50
100
150
200
250
300
300
50
A(i, j) = 1 if k=1:10((xi+k, yj+k)) > 7
Sequence 2
50
Sequence 1
100
150
200
250
300
50
150
100
Sequence 1
200
250
100
150
200
300
A dot matrix does not do any alignment
(global or local).
It helps to detect strongly conserved regions.
250
300
A(i, j) = 1 if k=1:20((xi+k, yj+k)) >
Sequence 2
50
50
Sequence 1
100
150
200
250
300
Seq1
Seq2
100
150
200
250
300
Today
• How to model gaps more accurately?
• Statistics of alignments
– Where does (xi, yj) come from?
– Are two aligned sequences actually related? –
not today
What’s a better alignment?
GACGCCGAACG
|||||
|||
GACGC---ACG
GACGCCGAACG
|||| | | ||
GACG-C-A-CG
Score = 8 x m – 3 x d
Score = 8 x m – 3 x d
However, gaps usually occur in bunches.
- During evolution, chunks of DNA may be lost entirely
- Aligning genomic sequences vs. cDNAs (reverse
complimentary to mRNAs)
Model gaps more accurately
• Current model:
– Gap of length n incurs penalty nd

n
• General:
– Convex function
– E.g. (n) = c * sqrt (n)

n
General gap dynamic
programming
Initialization:
same
Iteration:
F(i-1, j-1) + s(xi, yj)
F(i, j) = max maxk=0…i-1F(k,j) – (i-k)
maxk=0…j-1F(i,k) – (j-k)
Termination:
same
Running Time: O((M+N)MN)
(cubic)
Space: O(NM) (linear-space algorithm not applicable)
Compromise: affine gaps
(n) = d + (n – 1)e
|
|
gap
gap
open
extension
Match: 2
Gap open: -5
Gap extension: -1
(n)
d
GACGCCGAACG
|||||
|||
GACGC---ACG
GACGCCGAACG
|||| | | ||
GACG-C-A-CG
8x2-5-2 = 9
8x2-3x5 = 1
We want to find the optimal alignment with affine gap penalty in
• O(MN) time
• O(MN) or better O(M+N) memory
e
Allowing affine gap penalties
• Still three cases
– xi aligned with yj
– Xi aligns to a gap
• Are we continuing a gap in x? (if no, start is more expensive)
– Yj aligns to a gap
• Are we continuing a gap in y? (if no, start is more expensive)
• We can use a finite state machine to represent the three
cases as three states
– The machine has two heads, reading the chars on the two
strings separately
– At every step, each head reads 0 or 1 char from each sequence
– Depending on what it reads, goes to a different state, and
produces different scores
Finite State Machine
Input
Output
?/?
?/?
Ix
?/?
?/?
F
State
?/?
Iy
?/?
?/?
F: have just read 1 char from each seq (xi aligned to yj )
Ix: have read 0 char from x. (yj aligned to a gap)
Iy: have read 0 char from y (xi aligned to a gap)
Input
Output
(xi,yj) / 
(xi,yj) / 
Ix
(-, yj) / d
F
Start state
(-, yj) / e
(xi,-) / d
Iy
(xi,-) / e
(xi,yj) / 
Current state
Input
Output
Next state
F
(xi,yj)

F
F
(-,yj)
(xi,-)
(-,yj)
…
d
d
e
…
Ix
F
Ix
…
Iy
Ix
…
(xi,yj) / 
(xi,yj) /

start
state
(-, yj) / e
Ix
(-, yj) / d
F
(xi,-) / d
(xi,yj) / 
F-F-F-F
Iy
(xi,-) / e
F-Iy-F-F-Ix
AAC
AAC
AAC-
ACT
|||
||
ACT
-ACT
F-F-Iy-F-Ix
AAC-
| |
A-CT
Given a pair of sequences, an alignment (not necessarily optimal)
corresponds to a state path in the FSM.
Optimal alignment: find a state path to read the two sequences such that
the total output score is the highest
Dynamic programming
• We encode this information in three different
matrices
• For each element (i,j) we use three variables
– F(i,j): best alignment (score) of x1..xi & y1..yj if xi aligns
to yj
– Ix(i,j): best alignment of x1..xi & y1..yj if yj aligns to gap
– Iy(i,j): best alignment of x1..xi & y1..yj if xi aligns to gap
xi
xi
yj
F(i, j)
xi
yj
Ix(i, j)
yj
Iy(i, j)
(-, yj)/e
(xi,yj) /
(xi,yj) /
Ix
(-, yj) /d
F
(xi,-) /d
Iy
(xi,yj) /
xi
yj
(xi,-)/e
F(i-1, j-1) + (xi, yj)
F(i, j) = max Ix(i-1, j-1) + (xi, yj)
Iy(i-1, j-1) + (xi, yj)
(-, yj)/e
(xi,yj) /
(xi,yj) /
Ix
(-, yj) /d
F
(xi,-) /d
Iy
(xi,yj) /
(xi,-)/e
F(i, j-1) + d
xi
Ix(i, j) = max
yj
Ix(i, j)
Ix(i, j-1) + e
(-, yj)/e
(xi,yj) /
(xi,yj) /
Ix
(-, yj) /d
F
(xi,-) /d
Iy
(xi,yj) /
(xi,-)/e
F(i-1, j) + d
xi
Iy(i, j) = max
yj
Iy(i, j)
Iy(i-1, j) + e
F(i – 1, j – 1)
F(i, j) = (xi, yj) + max Ix(i – 1, j – 1)
Iy(i – 1, j – 1)
Ix(i, j) = max
Iy(i, j) = max
Continuing alignment
Closing gaps in x
Closing gaps in y
F(i, j – 1) + d
Opening a gap in x
Ix(i, j – 1) + e
Gap extension in x
F(i – 1, j) + d
Opening a gap in y
Iy(i – 1, j) + e
Gap extension in y
Data dependency
F
i
Ix
Iy
j
i-1
i-1
j-1
j-1
Data dependency
F
i
Ix
Iy
j
i
i
j
j
Data dependency
• If we stack all three matrices
– No cyclic dependency
– Therefore, we can fill in all three matrices in order
Algorithm
• for i = 1:m
– for j = 1:n
• Fill in F(i, j), Ix(i, j), Iy(i, j)
– end
end
• F(M, N) = max (F(M, N), Ix(M, N), Iy(M, N))
• Time: O(MN)
• Space: O(MN) or O(N) when combine with the
linear-space algorithm
Exercise
•
•
•
•
•
•
x = GCAC
y = GCC
m=2
s = -2
d = -5
e = -1
y= G
x= 0
-
C
-
C
-
x=
y= G
C
C
-
-
-
G -
G -5
C -
C -6
A -
A -7
C -
C -8
F: aligned on both
y= G
x = -
-5
C
-6
A
C
m=2
s = -2
d = -5
e = -1
Iy: Insertion on y
C
-7
F(i-1, j-1)
-
Iy(i-1, j-1)
(xi, yj)
G -
C
-
Iy(i-1,j)
F(i-1,j)
e
Ix(i-1, j-1)
d
F(i, j)
-
F(i,j-1)
-
Ix(i,j)
Ix(i,j-1)
Ix: Insertion on x
Iy(i,j)
d
e
y= G
x= 0
-
G -
C
-
C
-
2
x=
y= G
C
C
-
-
-
-
G -5
C -
C -6
A -
A -7
C -
C -8
F
y= G
x = -
-5
C
A
C
Iy
C
-6
C
-7
F(i-1, j-1)
G -
m=2
s = -2
d = -5
e = -1
Iy(i-1, j-1)
(xi, yj) = 2
-
Ix(i-1, j-1)
-
F(i, j)
-
Ix
y= G
x= 0
G -
C
-
-
2
-7
C
-
x=
y= G
C
C
-
-
-
-
G -5
C -
C -6
A -
A -7
C -
C -8
F
y= G
x = -
-5
C
A
C
Iy
C
-6
C
-7
F(i-1, j-1)
G -
m=2
s = -2
d = -5
e = -1
Iy(i-1, j-1)
(xi, yj) = -2
-
Ix(i-1, j-1)
-
F(i, j)
-
Ix
y= G
x= 0
G -
C
C
-
-
-
2
-7
-8
x=
y= G
C
C
-
-
-
-
G -5
C -
C -6
A -
A -7
C -
C -8
F
y= G
x = -
-5
C
A
C
Iy
C
-6
C
-7
F(i-1, j-1)
G -
m=2
s = -2
d = -5
e = -1
Iy(i-1, j-1)
(xi, yj) = -2
-
Ix(i-1, j-1)
-
F(i, j)
-
Ix
y= G
x= 0
G -
C
C
-
-
-
2
-7
-8
y= G
C
C
-
-
-
x=
G -5
C -
C -6
A -
A -7
C -
C -8
F
y= G
x=
G -
Iy
C
-5
-6
-
-3
C
-7
F(i,j-1)
C -
d = -5
Ix(i,j)
A -
Ix(i,j-1)
C -
Ix
e = -1
m=2
s = -2
d = -5
e = -1
y= G
x= 0
G -
C
C
-
-
-
2
-7
-8
y= G
C
C
-
-
-
x=
G -5
C -
C -6
A -
A -7
C -
C -8
F
y= G
x=
G -
Iy
C
C
-5
-6
-7
-
-3
-4
F(i,j-1)
C -
d = -5
Ix(i,j)
A -
Ix(i,j-1)
C -
Ix
e = -1
m=2
s = -2
d = -5
e = -1
y= G
x= 0
G -
C
C
-
-
-
2
-7
-8
y= G
C
C
-
-
-
-
-
-
x=
G -5
C -
C -6
A -
A -7
C -
C -8
F
y= G
x=
G -
-5
-
Iy
C
-6
-3
C
-7
-4
Iy(i-1,j)
F(i-1,j)
e=-1
C -
d=-5
A -
Iy(i,j)
C -
Ix
m=2
s = -2
d = -5
e = -1
y= G
x= 0
C
C
-
-
-
G -
2
-7
-8
C -
-7
y= G
C
C
-
-
-
-
-
-
x=
G -5
C -6
A -
A -7
C -
C -8
F
y= G
x=
G -
-5
-
m=2
s = -2
d = -5
e = -1
Iy
C
-6
-3
C
-7
-4
F(i-1, j-1)
Iy(i-1, j-1)
(xi, yj) = -2
C -
Ix(i-1, j-1)
A -
F(i, j)
C -
Ix
y= G
x= 0
C
C
-
-
-
G -
2
-7
-8
C -
-7
4
y= G
C
C
-
-
-
-
-
-
x=
G -5
C -6
A -
A -7
C -
C -8
F
y= G
x=
G -
-5
-
m=2
s = -2
d = -5
e = -1
Iy
C
-6
-3
C
-7
-4
F(i-1, j-1)
Iy(i-1, j-1)
(xi, yj) = 2
C -
Ix(i-1, j-1)
A -
F(i, j)
C -
Ix
y= G
x= 0
C
C
y= G
C
C
-
-
-
-
-
-
-
-
-
G -
2
-7
-8
G -5
C -
-7
4
-1
C -6
x=
A -
A -7
C -
C -8
F
y= G
x=
G -
-5
Iy
C
-6
-
m=2
s = -2
d = -5
e = -1
-3
C
-7
-4
F(i-1, j-1)
Iy(i-1, j-1)
(xi, yj) = 2
C -
Ix(i-1, j-1)
A -
F(i, j)
C -
Ix
y= G
x= 0
C
C
y= G
C
C
-
-
-
-
-
-
-
-
-
G -
2
-7
-8
G -5
C -
-7
4
-1
C -6
x=
A -
A -7
C -
C -8
F
y= G
x=
Iy
C
C
-5
-6
-7
G -
-
-3
-4
C -
-
-12 -1
F(i,j-1)
d = -5
Ix(i,j)
A -
Ix(i,j-1)
C -
Ix
e = -1
m=2
s = -2
d = -5
e = -1
y= G
x= 0
C
C
y= G
C
C
-
-
-
-
-
-
-
-
G -
2
-7
-8
G -5
-
C -
-7
4
-1
C -6
-3
x=
A -
A -7
C -
C -8
F
y= G
x=
-5
Iy
C
-6
C
-7
G -
-
-3
C -
-
-12 -1
-4
Iy(i-1,j)
F(i-1,j)
e=-1
d=-5
A -
Iy(i,j)
C -
Ix
m=2
s = -2
d = -5
e = -1
y= G
x= 0
C
C
y= G
C
C
-
-
-
-
-
-
-
G -
2
-7
-8
G -5
-
-
C -
-7
4
-1
C -6
-3
-12 -13
x=
A -
A -7
C -
C -8
F
y= G
x=
G -
Iy
C
C
-5
-6
-7
-
-3
-4
F(i-1, j-1)
Iy(i-1, j-1)
(xi, yj)
-
Iy(i-1,j)
F(i-1,j)
e
Ix(i-1, j-1)
C -
m=2
s = -2
d = -5
e = -1
d
-12 -1
F(i, j)
F(i,j-1)
A -
Iy(i,j)
d
Ix(i,j)
C -
Ix(i,j-1)
Ix
e
y= G
x= 0
C
C
y= G
C
C
-
-
-
-
-
-
-
G -
2
-7
-8
G -5
-
-
C -
-7
4
-5
C -6
-3
-12 -13
A -
-8
-5
2
A -7
C -
x=
C -8
F
y= G
x=
G -
m=2
s = -2
d = -5
e = -1
Iy
C
C
-5
-6
-7
-
-3
-4
F(i-1, j-1)
Iy(i-1, j-1)
(xi, yj)
F(i-1,j)
e
Ix(i-1, j-1)
C -
-
A -
-
Iy(i-1,j)
d
-12 -1
F(i, j)
-13 -10
F(i,j-1)
Iy(i,j)
d
Ix(i,j)
C -
Ix(i,j-1)
Ix
e
y= G
x= 0
C
C
y= G
C
C
-
-
-
-
-
-
-
G -
2
-7
-8
G -5
-
-
C -
-7
4
-1
C -6
-3
-12 -13
A -
-8
-5
2
A -7
-8
-1
C -
x=
C -8
F
y= G
x=
-5
Iy
C
-6
C
-7
G -
-
-3
C -
-
-12 -1
A -
-
-4
Iy(i-1,j)
F(i-1,j)
e=-1
d=-5
-13 -10
Iy(i,j)
C -
Ix
m=2
s = -2
d = -5
e = -1
y= G
x= 0
C
C
y= G
C
C
-
-
-
-
-
-
-
G -
2
-7
-8
G -5
-
-
C -
-7
4
-1
C -6
-3
-12 -13
A -
-8
-5
2
A -7
-8
-1
C -
x=
C -8
F
y= G
x=
-6
-5
Iy
C
-6
C
-7
G -
-
-3
C -
-
-12 -1
A -
-
-4
Iy(i-1,j)
F(i-1,j)
e=-1
d=-5
-13 -10
Iy(i,j)
C -
Ix
m=2
s = -2
d = -5
e = -1
y= G
x= 0
C
C
y= G
C
C
-
-
-
-
-
-
-
G -
2
-7
-8
G -5
-
-
C -
-7
4
-1
C -6
-3
-12 -13
A -
-8
-5
2
A -7
-8
-1
-6
C -
-9
-6
1
C -8
-13 -2
-3
x=
F
y= G
x=
G -
m=2
s = -2
d = -5
e = -1
Iy
C
C
-5
-6
-7
-
-3
-4
F(i-1, j-1)
Iy(i-1, j-1)
(xi, yj)
F(i-1,j)
e
Ix(i-1, j-1)
C -
-
A -
-
-13 -10
C -
-
-14 -11
Iy(i-1,j)
d
-12 -1
F(i, j)
F(i,j-1)
Ix(i,j)
Ix(i,j-1)
Ix
Iy(i,j)
d
e
y= G
x= 0
C
C
y= G
C
C
-
-
-
-
-
-
-
G -
2
-7
-8
G -5
-
-
C -
-7
4
-1
C -6
-3
-12 -13
A -
-8
-5
2
A -7
-8
-1
-6
C -
-9
-6
1
C -8
-13 -2
-3
x=
F
y= G
x=
G -
m=2
s = -2
d = -5
e = -1
Iy
C
C
-5
-6
-7
-
-3
-4
F(i-1, j-1)
Iy(i-1, j-1)
(xi, yj)
F(i-1,j)
e
Ix(i-1, j-1)
C -
-
A -
-
-13 -10
C -
-
-14 -11
Iy(i-1,j)
d
-12 -1
F(i, j)
F(i,j-1)
Ix(i,j)
Ix(i,j-1)
Ix
Iy(i,j)
d
e
y= G
x= 0
C
C
y= G
C
C
-
-
-
-
-
-
-
G -
2
-7
-8
G -5
-
-
C -
-7
4
-1
C -6
-3
-12 -13
A -
-8
-5
2
x GCAC
A -7
-8
-1
-6
C -
-9
-6
1
|| |
C -8
-13 -2
-3
y GC-C
F
y= G
x=
x=
C
Iy
C
y= G
-5
-6
-7
G -
-
-3
-4
G
C -
-
-12 -1
C
A -
-
-13 -10
C -
-
-14 -11
Ix
m=2
s = -2
d = -5
e = -1
x=
A
C
C
C