Transcript Fasta, Blast, Probabilities

Fasta, Blast, Probabilities
.
Reminder
 Last
classes we discussed dynamic programming
algorithms for
 global alignment
 local alignment
 Multiple alignment
 All
of these assumed a scoring rule:
 : (  {})  (  {})  
that determines the quality of perfect matches,
substitutions, insertions, and deletions.
2
Alignment in Real Life
 One
of the major uses of alignments is to find
sequences in a “database.”
 The current protein database contains about 100
millions (i.e.,108) residues ! So searching a 1000
long target sequence requires to evaluate about
1011 matrix cells which will take about three hours in
the rate of 10 millions evaluations per second.
 Quite annoying when, say, one thousand target
sequences need to be searched because it will take
about four months to run.
3
Heuristic Fast Search
 Instead,
most searches rely on heuristic procedures
 These are not guaranteed to find the best match
 Sometimes, they will completely miss a high-scoring
match
We now describe the main ideas used by the best
known of these heuristic procedures.
4
Basic Intuition
 Almost
all heuristic search procedures are based
on the observation that real-life matches often
contain long strings with gap-less matches.
 These
heuristic try to find significant gap-less
matches and then extend them.
5
Banded DP
that we have two strings s[1..n] and
t[1..m] such that nm
 If the optimal alignment of s and t has few gaps,
then path of the alignment will be close to diagonal
 Suppose
s
t
6
Banded DP
 To
find such a path, it suffices to search in a
diagonal region of the matrix.
 If the diagonal band has width k, then the dynamic
programming step takes O(kn).
 Much faster than O(n2) of standard DP.
s
V[i,i+k/2]
V[i, i+k/2+1]
Out of range
V[i+1, i+k/2 +1]
Note that for diagonals i-j = constant.
t
k
7
Banded DP for local alignment
Problem: Where is the banded diagonal ? It need
not be the main diagonal when looking for a good
local alignment.
s
How do we select which
subsequences to align
using banded DP?
k
t
We heuristically find potential
diagonals and evaluate them
using Banded DP.
This is the main idea of FASTA.
8
Finding Potential Diagonals
Suppose that we have a relatively long gap-less match
AGCGCCATGGATTGAGCGA
TGCGACATTGATCGACCTA
Can we find “clues” that will let us find it quickly?
 Each such sequence defines a potential diagonal (which is
then evaluated using Banded DP.

9
Signature of a Match
Assumption: good matches contain several “patches” of
perfect matches
AGCGCCATGGATTGAGCTA
TGCGACATTGATCGACCTA
s
Since this is a gap-less alignment,
all perfect match regions
should be on one diagonal
t
10
FASTA-finding ungapped matches
Input: strings s and t, and a parameter ktup
 Find all pairs (i,j) such that s[i..i+ktup]=t[j..j+ktup]
 Locate sets of pairs that are on the same diagonal
 By sorting according to the difference i-j
 Compute the score for the diagonal that contains all these
pairs
s
t
11
FASTA-finding ungapped matches
Input: strings s and t, and a parameter ktup
 Find all pairs (i,j) such that s[i..i+ktup]=t[j..j+ktup]
 Step one: prepare an index of the database and the
query sequence such that given a sequence of
length ktup, one gets the list of positions. (Linear
time).
 Step two: for each ktup from the query add 1 in the
diagonal (i-j) in which it appears. Then find
contiguous (possibly with mismatch) ktup in s
diagonals.
t
12
FASTA- using banded DP
Step 3:
 Select the ten high scoring contiguous
segments
 Try and score all combinations of
these ten segments in order to
constitute a pass into the matrix
Step 4:
 Run banded DP on the region
containing the best scoring pass (say
with width 12).
s
2
Hence, the algorithm may combine
some diagonals into gapped
t
matches (in the example below
combine diagonals 2 and 3).
3
1
13
FASTA- practical choices
Most applications of FASTA use very small ktup (1-2 for
proteins, and 4-6 for DNA).
Higher values are faster, yielding less diagonal to search
around, but increase the chance to miss the optimal local
alignment.
s
t
Some implementation choices /tricks
have not been explicated herein.
14
FASTA-summary
Input: strings s and t, and a parameter ktup = 1,2,4,5, or 6
depending on the application.
Output: A highly scored local alignment
1.
2.
3.
Find pairs of matching substrings s[i..i+ktup]=t[j..j+ktup]
Extend to ungapped diagonals
Extend to gapped matches using banded DP
15
BLAST Overview
(Basic Local Alignment Search Tool)
Input: strings s and t, and a parameter T = threshold value
Output: A highly scored local alignment
Definition: Two strings s and t of length k are a high scoring
pair (HSP) if d(s,t) > T (usually consider un-gapped
alignments only).
1.
2.
3.
Find high scoring pairs of substrings such that d(s,t) > T

These words serve as seeds for finding longer matches
Extend to ungapped diagonals (as in FASTA)
Extend to gapped matches
16
BLAST Overview (cont.)
Step 1: Find high scoring pairs of substrings such that
d(s,t) > T (The seeds):

Find all strings of length k which score at
least T with substrings of s in a gapless
alignment (k = 4 for proteins, 11 for DNA)
(note: possibly, not all k-words must be tested, e.g. when
such a word scores less than T with itself).

Find in t all exact matches with each of the
above strings.
17
Extending Potential Matches
Once a seed is found, BLAST attempts to find
a local alignment that extends the seed.
Seeds on the same diagonal are combined
(as in FASTA), then extended as far as
possible in a greedy manner without gap.
During the extension phase, the
search stops when the score passes
below some lower bound computed by
BLAST (to save time).
For the best ungap alignment do a
banded SW an assign a probabilistic
score.
s
t
18
19
20
21
22
23
24
25
Why use probability to define and/or
interpret a scoring function ?
• Similarity is probabilistic in nature because
biological changes like mutation, recombination, and
selection, are not deterministic.
• We could answer questions such as:
• How probable two sequences are similar?
• Is the similarity found significant or random?
• How to change a similarity score when, say,
mutation rate of a specific area on the
chromosome becomes known ?
26
A Probabilistic Model
 For
now, we will focus on alignment without indels.
 For now, we assume each position (nucleotide
/amino-acid) is independent of other positions.
 We consider two options:
M: the sequences are Matched (related)
R: the sequences are Random (unrelated)
27
Unrelated Sequences
 Our
random model of unrelated sequences is simple
 Each position is sampled independently from a
distribution over the alphabet 
 We assume there is a distribution q() that
describes the probability of letters in such
positions.
 Then:
P( s[1..n], t[1..n] | R)   q( s[i]) q(t[i])
i
28
Related Sequences
 We
assume that each pair of aligned positions
(s[i],t[i]) evolved from a common ancestor
 Let p(a,b) be a distribution over pairs of letters.
 p(a,b) is the probability that some ancestral letter
evolved into this particular pair of letters.
P( s[1..n], t[1..n] | M )   p( s[i], t[i])
i
29
Odd-Ratio Test for Alignment
p( s[i ], t[i ])

P ( s, t | M )
p ( s[i ], t[i ])
Q


P( s, t | R)  q ( s[i ]) q (t[i ])
q( s[i ]) q (t[i ])
i
i
i
If Q > 1, then the two strings s and t are more likely to
be related (M) than unrelated (R).
If Q < 1, then the two strings s and t are more likely to
be unrelated (R) than related (M).
30
Log Odd-Ratio Test for Alignment
Taking logarithm of Q yields
Score(s[i],t[i])
P ( s, t | M )
p ( s[i ], t[i ])
p ( s[i ], t[i ])
log
 log 
  log
P ( s, t | R )
q ( s[i ]) q (t[i ])
i
i q ( s[i ]) q (t[i ])
If log Q > 0, then s and t are more likely to be related.
If log Q < 0, then they are more likely to be unrelated.
How can we relate this quantity to a score function ?
31
Probabilistic Interpretation of Scores
 We
define the scoring function via
p (a , b )
 (a , b )  log
q (a )q (b )
 Then,
the score of an alignment is the log-ratio
between the two models:

Score > 0

Score < 0


Model is more likely
Random is more likely
32
Estimating Probabilities
we are given a long string s[1..n] of
letters from 
 We want to estimate the distribution q(·) that
generated the sequence
 How should we go about this?
 Suppose
We build on the theory of parameter estimation in
statistics using either maximum likelihood
estimation or the Bayesian approach .
33
Estimating q()
 Suppose
from 

we are given a long string s[1..n] of letters
s can be the concatenation of all sequences in our
database
 We
want to estimate the distribution q()
n
Na
L
(
q
|
s
)

q
(
s
[
i
])

q
(
a
)
Likelihood function:


i 1
 That
a
is, q is defined per letter
34
Estimating q() (cont.)
Likelihood function: L(q | s ) 
How do we define q?
n
 q(s[i])   q(a)
i 1
a
ML parameters
MAP parameters
(Maximum Likelihood)
(Maximum A posteriori Probability)
Na
q(a) 
n
Na
Na 1
q(a) 
n |  |
35
Estimating p(·,·)
Intuition:
 Find pair of aligned sequences s[1..n], t[1..n],
 Estimate probability of pairs:
p (a , b ) 
Na ,b
n
Number of times a is
aligned with b in (s,t)
 Again,
s and t can be the concatenation of many
aligned pairs from the database
36
Problems in Estimating p(·,·)
 How
do we find pairs of aligned sequences?
 How far is the ancestor ?
 earlier divergence  low sequence similarity
 later divergence  high sequence similarity
 Does one letter mutate to the other or are they both
mutations of a common ancestor having yet
another residue/nucleotide acid ?
37
47
BLOSUM 62
50