Transcript Lecture_6-2

An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Sequence Alignment
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Outline
•
•
•
•
Global Alignment
Scoring Matrices
Local Alignment
Alignment with Affine Gap Penalties
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
From LCS to Alignment: Change Scoring
• The Longest Common Subsequence (LCS)
problem—the simplest form of sequence
alignment – allows only insertions and deletions (no
mismatches).
• In the LCS Problem, we scored 1 for matches and 0
for indels
• Consider penalizing indels and mismatches with
negative scores
• Simplest scoring schema:
+1 : match premium
-μ : mismatch penalty
-σ : indel penalty
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Simple Scoring
• When mismatches are penalized by –μ,
indels are penalized by –σ,
and matches are rewarded with +1,
the resulting score is:
#matches – μ(#mismatches) – σ (#indels)
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
The Global Alignment Problem
Find the best alignment between two strings under a given scoring
schema
Input : Strings v and w and a scoring schema
Output : Alignment of maximum score
→ = -б
= 1 if match
= -µ if mismatch
si,j = max
si-1,j-1 +1 if vi = wj
s i-1,j-1 -µ if vi ≠ wj
s i-1,j - σ
s i,j-1 - σ
m : mismatch penalty
σ
: indel penalty
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Scoring Matrices
To generalize scoring, consider a (4+1) x(4+1) scoring
matrix δ.
In the case of an amino acid sequence alignment, the
scoring matrix would be a (20+1)x(20+1) size. The
addition of 1 is to include the score for comparison
of a gap character “-”.
This will simplify the algorithm as follows:
si-1,j-1 + δ (vi, wj)
si,j = max
s i-1,j + δ (vi, -)
s i,j-1 + δ (-, wj)
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Measuring Similarity
• Measuring the extent of similarity between
two sequences
• Based on percent sequence identity
• Based on conservation
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Percent Sequence Identity
• The extent to which two nucleotide or amino
acid sequences are invariant
AC C TG A G – AG
AC G TG – G C AG
mismatch
indel
70% identical
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Making a Scoring Matrix
• Scoring matrices are created based on
biological evidence.
• Alignments can be thought of as two
sequences that differ due to mutations.
• Some of these mutations have little effect on
the protein’s function, therefore some
penalties, δ(vi , wj), will be less harsh than
others.
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Scoring Matrix: Example
A
R
N
K
A
5
-2
-1
-1
R
-
7
-1
3
N
-
-
7
0
K
-
-
-
6
• Notice that although
R and K are different
amino acids, they
have a positive score.
• Why? They are both
positively charged
amino acids will not
greatly change
function of protein.
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Scoring matrices
• Amino acid substitution matrices
• PAM
• BLOSUM
• DNA substitution matrices
• DNA is less conserved than protein
sequences
• Less effective to compare coding regions at
nucleotide level
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
PAM
• Point Accepted Mutation (Dayhoff et al.)
• 1 PAM = PAM1 = 1% average change of all amino
acid positions
• After 100 PAMs of evolution, not every residue will
have changed
• some residues may have mutated several
times
• some residues may have returned to their
original state
• some residues may not changed at all
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
PAMX
• PAMx = PAM1x
• PAM250 = PAM1250
• PAM250 is a widely used scoring matrix:
Ala
Arg
Asn
Asp
Cys
Gln
...
Trp
Tyr
Val
A
R
N
D
C
Q
Ala
A
13
3
4
5
2
3
Arg
R
6
17
4
4
1
5
Asn
N
9
4
6
8
1
5
Asp
D
9
3
7
11
1
6
Cys
C
5
2
2
1
52
1
Gln
Q
8
5
5
7
1
10
Glu
E
9
3
6
10
1
7
Gly
G
12
2
4
5
2
3
His
H
6
6
6
6
2
7
Ile
I
8
3
3
3
2
2
W
Y
V
0
1
7
2
1
4
0
2
4
0
1
4
0
3
4
0
1
4
0
1
4
0
1
4
1
3
5
0
2
4
Leu
L
6
2
2
2
1
3
Lys ...
K ...
7 ...
9
5
5
1
5
1
2
15
0
1
10
A
R
N
D
C
Q
E
G
H
I
L
K
M
F
P
S
T
W
Y
V
2
-2 6
0 0 2
0 -1 2 4
-2 -4 -4 -5 12
0 1 1 2 -5 4
0 -1 1 3 -5 2 4
1 -3 0 1 -3 -1 0 5
-1 2 2 1 -3 3 1 -2 6
-1 -2 -2 -2 -2 -2 -2 -3 -2 5
-2 -3 -3 -4 -6 -2 -3 -4 -2 -2 6
-1 3 1 0 -5 1 0 -2 0 -2 -3 5
-1 0 -2 -3 -5 -1 -2 -3 -2 2 4 0 6
-3 -4 -3 -6 -4 -5 -5 -5 -2 1 2 -5 0 9
1 0 0 -1 -3 0 -1 0 0 -2 -3 -1 -2 -5 6
1 0 1 0 0 -1 0 1 -1 -1 -3 0 -2 -3 1 2
1 -1 0 0 -2 -1 0 0 -1 0 -2 0 -1 -3 0 1 3
-6 2 -4 -7 -8 -5 -7 -7 -3 -5 -2 -3 -4 0 -6 -2 -5 17
-3 -4 -2 -4 0 -4 -4 -5 0 -1 -1 -4 -2 7 -5 -3 -3 0 10
0 -2 -2 -2 -2 -2 -2 -1 -2 4 2 -2 2 -1 -1 -1 0 -6 -2 4
A R N D C Q E G H I
L K M F P S T W Y V
PAM250 log odds
scoring matrix
Why do we go from a mutation probability
matrix to a log odds matrix?
• We want a scoring matrix so that when we do a pairwise
alignment (or a BLAST search) we know what score to
assign to two aligned amino acid residues.
• Logarithms are easier to use for a scoring system. They
allow us to sum the scores of aligned residues (rather
than having to multiply them).
How do we go from a mutation probability
matrix to a log odds matrix?
• The cells in a log odds matrix consist of an “odds ratio”:
the probability that an alignment is authentic
the probability that the alignment was random
The score S for an alignment of residues a,b is given by:
S(a,b) = 10 log10 ( Mab / pb )
As an example, for tryptophan,
S( W, W ) = 10 log10 ( 0.55 / 0.01 ) = 17.4
What do the numbers mean
in a log odds matrix?
S( W, W ) = 10 log10 ( 0.55 / 0.010 ) = 17.4
A score of +17 for tryptophan means that this alignment
is 50 times more likely than a chance alignment of two
tryptophan residues.
S(W, W) = 17
Probability of replacement ( Mab / pb ) = x
Then
17 = 10 log10 x
1.7 = log10 x
101.7 = x = 50
What do the numbers mean
in a log odds matrix?
A score of +2 indicates that the amino acid replacement
occurs 1.6 times as frequently as expected by chance.
A score of 0 is neutral.
A score of –10 indicates that the correspondence of two
amino acids in an alignment that accurately represents
homology (evolutionary descent) is one tenth as frequent
as the chance alignment of these amino acids.
A
R
N
D
C
Q
E
G
H
I
L
K
M
F
P
S
T
W
Y
V
7
-10
9
-7
-9
9
-6
-17
-1
8
-10
-11
-17
-21
10
-7
-4
-7
-6
-20
9
-5
-15
-5
0
-20
-1
8
-4
-13
-6
-6
-13
-10
-7
7
-11
-4
-2
-7
-10
-2
-9
-13
10
-8
-8
-8
-11
-9
-11
-8
-17
-13
9
-9
-12
-10
-19
-21
-8
-13
-14
-9
-4
7
-10
-2
-4
-8
-20
-6
-7
-10
-10
-9
-11
7
-8
-7
-15
-17
-20
-7
-10
-12
-17
-3
-2
-4
12
-12
-12
-12
-21
-19
-19
-20
-12
-9
-5
-5
-20
-7
9
-4
-7
-9
-12
-11
-6
-9
-10
-7
-12
-10
-10
-11
-13
8
-3
-6
-2
-7
-6
-8
-7
-4
-9
-10
-12
-7
-8
-9
-4
7
-3
-10
-5
-8
-11
-9
-9
-10
-11
-5
-10
-6
-7
-12
-7
-2
8
-20
-5
-11
-21
-22
-19
-23
-21
-10
-20
-9
-18
-19
-7
-20
-8
-19
13
-11
-14
-7
-17
-7
-18
-11
-20
-6
-9
-10
-12
-17
-1
-20
-10
-9
-8
10
-5
-11
-12
-11
-9
-10
-10
-9
-9
-1
-5
-13
-4
-12
-9
-10
-6
-22
-10
R
N
D
Q
E
A
C
G
H
PAM10 log odds
scoring matrix
I
L
K
M
F
P
S
T
W Y
8
V
BLOSUM Outline
• Idea: use aligned ungapped regions of protein
families.These are assumed to have a common
ancestor. Similar ideas but better statistics and
modeling.
• Procedure:
– Cluster together sequences in a family whenever more than
L% identical residues are shared, for BLOSUM-L.
– Count number of substitutions across different clusters (in
the same family).
– Estimate frequencies using the counts.
• Practice: BlOSUM-50 and BLOSOM62 are wildly
used.
Considered the state of the art nowadays.
BLOSUM Matrices
BLOSUM matrices are based on local alignments.
BLOSUM stands for blocks substitution matrix.
BLOSUM62 is a matrix calculated from comparisons of
sequences with less than 62% identical sites.
BLOSUM Matrices
Percent amino acid identity
100
62
30
BLOSUM62
Percent amino acid identity
BLOSUM Matrices
100
100
100
62
62
62
30
30
30
BLOSUM80
BLOSUM62
BLOSUM30
BLOSUM Matrices
All BLOSUM matrices are based on observed alignments;
they are not extrapolated from comparisons of
closely related proteins.
The BLOCKS database contains thousands of groups of
multiple sequence alignments.
BLOSUM62 is the default matrix in BLAST 2.0.
Though it is tailored for comparisons of moderately distant
proteins, it performs well in detecting closer relationships.
A search for distant relatives may be more sensitive
with a different matrix.
A
R
N
D
C
Q
E
G
H
I
L
K
M
F
P
S
T
W
Y
V
4
-1 5
-2 0 6
-2 -2 1 6
0 -3 -3 -3 9
-1 1 0 0 -3 5
-1 0 0 2 -4 2 5
0 -2 0 -1 -3 -2 -2 6
-2 0 1 -1 -3 0 0 -2 8
-1 -3 -3 -3 -1 -3 -3 -4 -3 4
-1 -2 -3 -4 -1 -2 -3 -4 -3 2 4
-1 2 0 -1 -1 1 1 -2 -1 -3 -2 5
-1 -2 -2 -3 -1 0 -2 -3 -2 1 2 -1 5
-2 -3 -3 -3 -2 -3 -3 -3 -1 0 0 -3 0 6
-1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4 7
1 -1 1 0 -1 0 0 0 -1 -2 -2 0 -1 -2 -1 4
0 -1 0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1 1 5
-3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1 1 -4 -3 -2 11
-2 -2 -2 -3 -2 -1 -2 -3 2 -1 -1 -2 -1 3 -3 -2 -2
2
7
0 -3 -3 -3 -1 -2 -2 -3 -3 3 1 -2 1 -1 -2 -2 0 -3 -1 4
A R N D C Q E G H I
L K M F P S T W Y
V
Blosum62 scoring matrix
Rat versus
mouse RBP
Rat versus
bacterial
lipocalin
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local vs. Global Alignment
• The Global Alignment Problem tries to find
the longest path between vertices (0,0) and
(n,m) in the edit graph.
• The Local Alignment Problem tries to find the
longest path among paths between arbitrary
vertices (i,j) and (i’, j’) in the edit graph.
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local vs. Global Alignment
• The Global Alignment Problem tries to find the
longest path between vertices (0,0) and (n,m) in the
edit graph.
• The Local Alignment Problem tries to find the
longest path among paths between arbitrary
vertices (i,j) and (i’, j’) in the edit graph.
• In the edit graph with negatively-scored edges,
Local Alignmet may score higher than Global
Alignment
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local vs. Global Alignment (cont’d)
• Global Alignment
--T—-CC-C-AGT—-TATGT-CAGGGGACACG—A-GCATGCAGA-GAC
| || | || | | | |||
|| | | | | ||||
|
AATTGCCGCC-GTCGT-T-TTCAG----CA-GTTATG—T-CAGAT--C
• Local Alignment—better alignment to find
conserved segment
tccCAGTTATGTCAGgggacacgagcatgcagagac
||||||||||||
aattgccgccgtcgttttcagCAGTTATGTCAGatc
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local Alignment: Example
Local alignment
Global alignment
Compute a “mini”
Global Alignment to
get Local
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local Alignments: Why?
• Two genes in different species may be similar
over short conserved regions and dissimilar
over remaining regions.
• Example:
• Homeobox genes have a short region
called the homeodomain that is highly
conserved between species.
• A global alignment would not find the
homeodomain because it would try to align
the ENTIRE sequence
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
The Local Alignment Problem
• Goal: Find the best local alignment between
two strings
• Input : Strings v, w and scoring matrix δ
• Output : Alignment of substrings of v and w
whose alignment score is maximum among
all possible alignment of all possible
substrings
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
The Problem with this Problem
• Long run time O(n4):
- In the grid of size n x n there are ~n2
vertices (i,j) that may serve as a source.
- For each such vertex computing alignments
from (i,j) to (i’,j’) takes O(n2) time.
• This can be remedied by giving free rides
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local Alignment: Example
Local alignment
Global alignment
Compute a “mini”
Global Alignment to
get Local
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local Alignment: Example
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local Alignment: Example
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local Alignment: Example
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local Alignment: Example
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local Alignment: Example
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local Alignment: Running Time
• Long run time O(n4):
- In the grid of size n x n
there are ~n2 vertices (i,j)
that may serve as a
source.
- For each such vertex
computing alignments
from (i,j) to (i’,j’) takes
O(n2) time.
• This can be remedied by
giving free rides
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Local Alignment: Free Rides
Yeah, a free ride!
Vertex (0,0)
The dashed edges represent the free rides from
(0,0) to every other node.
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
The Local Alignment Recurrence
• The largest value of si,j over the whole edit
graph is the score of the best local alignment.
• The recurrence:
0
si,j = max
si-1,j-1 + δ (vi, wj)
s i-1,j + δ (vi, -)
s i,j-1 + δ (-, wj)
Notice there is only
this change from the
original recurrence of
a Global Alignment
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
The Local Alignment Recurrence
• The largest value of si,j over the whole edit
graph is the score of the best local alignment.
• The recurrence:
0
si,j = max
si-1,j-1 + δ (vi, wj)
s i-1,j + δ (vi, -)
s i,j-1 + δ (-, wj)
Power of ZERO: there is
only this change from the
original recurrence of a
Global Alignment - since
there is only one “free ride”
edge entering into every
vertex
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Scoring Indels: Naive Approach
• A fixed penalty σ is given to every indel:
• -σ for 1 indel,
• -2σ for 2 consecutive indels
• -3σ for 3 consecutive indels, etc.
Can be too severe penalty for a series of
100 consecutive indels
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Affine Gap Penalties
• In nature, a series of k indels often come as a
single event rather than a series of k single
nucleotide events:
This is more
likely.
Normal scoring would
give the same score This is less
for both alignments
likely.
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Accounting for Gaps
• Gaps- contiguous sequence of spaces in one of the
rows
• Score for a gap of length x is:
-(ρ + σx)
where ρ >0 is the penalty for introducing a gap:
gap opening penalty
ρ will be large relative to σ:
gap extension penalty
because you do not want to add too much of a
penalty for extending the gap.
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Affine Gap Penalties
• Gap penalties:
• -ρ-σ when there is 1 indel
• -ρ-2σ when there are 2 indels
• -ρ-3σ when there are 3 indels, etc.
• -ρ- x·σ (-gap opening - x gap extensions)
• Somehow reduced penalties (as compared to
naïve scoring) are given to runs of horizontal
and vertical edges
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Affine Gap Penalties and Edit Graph
To reflect affine gap
penalties we have to
add “long” horizontal
and vertical edges to
the edit graph. Each
such edge of length x
should have weight
- - x *
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Adding “Affine Penalty” Edges to the Edit Graph
There are many such edges!
Adding them to the graph
increases the running time
of the alignment algorithm
by a factor of n (where n is
the number of vertices)
So the complexity increases
from O(n2) to O(n3)
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Manhattan in 3 Layers
ρ
δ
δ
σ
δ
ρ
δ
δ
σ
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Affine Gap Penalties and 3 Layer Manhattan Grid
• The three recurrences for the scoring
algorithm creates a 3-layered graph.
• The top level creates/extends gaps in the
sequence w.
• The bottom level creates/extends gaps in
sequence v.
• The middle level extends matches and
mismatches.
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Switching between 3 Layers
• Levels:
• The main level is for diagonal edges
• The lower level is for horizontal edges
• The upper level is for vertical edges
• A jumping penalty is assigned to moving from the
main level to either the upper level or the lower level
(-- )
• There is a gap extension penalty for each
continuation on a level other than the main level (-)
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
The 3 Grids
-σ
-(ρ + σ)
+δ(vi,wj)
+0
+0
-σ
Away from mid-level
Toward mid-level
-(ρ + σ)
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
The 3-leveled Manhattan Grid
An Introduction to Bioinformatics Algorithms
www.bioalgorithms.info
Affine Gap Penalty Recurrences
si,j =
max
s i-1,j - σ
s i-1,j –(ρ+σ)
Continue Gap in w (deletion)
Start Gap in w (deletion): from middle
si,j =
max
s i,j-1 - σ
s i,j-1 –(ρ+σ)
Continue Gap in v (insertion)
si,j =
max
si-1,j-1 + δ (vi, wj) Match or Mismatch
s i,j
End deletion: from top
s i,j
End insertion: from bottom
Start Gap in v (insertion):from middle