- 93 Views
- Uploaded on
- Presentation posted in: General

Sequence Alignment I Lecture #2

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Sequence Alignment ILecture #2

Background Readings: The second chapter (pages 12-45) in the text book, Biological Sequence Analysis, Durbin et al., 2001.

This class has been edited from Nir Friedman’s lecture which is available at www.cs.huji.ac.il/~nir. Changes made by Dan Geiger. Thanks to Shlomo Moran for several improvements.

.

Much of bioinformatics involves sequences

- DNA sequences
- RNA sequences
- Protein sequences
We can think of these sequences as strings of letters

- DNA & RNA: alphabet of 4 letters
- Protein: alphabet of 20 letters

g1

g2

- Finding similarity between sequences is important for many biological questions
For example:

- Find similar proteins
- Allows to predict function & structure

- Locate similar subsequences in DNA
- Allows to identify (e.g) regulatory elements

- Locate DNA sequences that might overlap
- Helps in sequence assembly

Input: two sequences over the same alphabet

Output: an alignment of the two sequences

Example:

- GCGCATGGATTGAGCGA
- TGCGCCATTGATGACCA
A possible alignment:

-GCGC-ATGGATTGAGCGA

TGCGCCATTGAT-GACC-A

-GCGC-ATGGATTGAGCGA

TGCGCCATTGAT-GACC-A

Three elements:

- Perfect matches
- Mismatches
- Insertions & deletions (indel)

There are many possible alignments

For example, compare:

-GCGC-ATGGATTGAGCGA

TGCGCCATTGAT-GACC-A

to

------GCGCATGGATTGAGCGA

TGCGCC----ATTGATGACCA--

Which one is better?

Intuition:

- Similar sequences evolved from a common ancestor
- Evolution changed the sequences from this ancestral sequence by mutations:
- Replacements: one letter replaced by another
- Deletion: deletion of a letter
- Insertion: insertion of a letter

- Scoring of sequence similarity should examine how many and which operations took place

Score each position independently:

- Match: +1
- Mismatch: -1
- Indel -2
Score of an alignment is sum of position scores

Example:

-GCGC-ATGGATTGAGCGA

TGCGCCATTGAT-GACC-A

Score: (+1x13) + (-1x2) + (-2x4) = 3

------GCGCATGGATTGAGCGA

TGCGCC----ATTGATGACCA--

Score: (+1x5) + (-1x6) + (-2x11) = -23

- The choice of +1,-1, and -2 scores is quite arbitrary
- Depending on the context, some changes are moreplausible than others
- Exchange of an amino-acid by one with similar properties (size, charge, etc.)
vs.

- Exchange of an amino-acid by one with opposite properties

- Exchange of an amino-acid by one with similar properties (size, charge, etc.)
- Probabilistic interpretation: (e.g.) How likely is one alignment versus another ?

- We define a scoring function by specifying a function
- (x,y) is the score of replacing x by y
- (x,-) is the score of deleting x
- (-,x) is the score of inserting x

- The score of an alignment is the sum of position scores

- The optimal (maximal) score between two sequences is the maximal score of all alignments of these sequences, namely,
- Computing the maximal score or actually finding an alignment that yields the maximal score are closely related tasks with similar algorithms.
- We now address these problems.

- How can we compute the optimal score ?
- If |s| = n and |t| = m, there are many alignments !
- Exercise : Show that the number of legal alignments denoted by A(m,n), where n > m, satisfies (legal means not to a “-” versus another “-”) :

- The additive form of the score allows us to perform dynamic programming to compute optimal score efficiently.

- Suppose we have two sequences:s[1..n+1] and t[1..m+1]
The best alignment must be one of three cases:

1. Last match is (s[n+1],t[m +1] )

2. Last match is (s[n +1],-)

3. Last match is (-, t[m +1] )

- Suppose we have two sequences:s[1..n+1] and t[1..m+1]
The best alignment must be one of three cases:

1. Last match is (s[n+1],t[m +1] )

2. Last match is (s[n +1],-)

3. Last match is (-, t[m +1] )

- Suppose we have two sequences:s[1..n+1] and t[1..m+1]
The best alignment must be one of three cases:

1. Last match is (s[n+1],t[m +1] )

2. Last match is (s[n +1],-)

3. Last match is (-, t[m +1] )

Define the notation:

- Using our recursive argument, we get the following recurrence for V:

T

S

- Of course, we also need to handle the base cases in the recursion:

AA

- -

versus

We fill the matrix using the recurrence rule:

T

S

We continue to fill the matrix using the recurrence rule

T

S

+1

-2 -A

A-

-2 (A- versus -A)

versus

T

S

T

S

Conclusion: d(AAAC,AGC) = -1

T

S

- To reconstruct the best alignment, we record which case(s) in the recursive rule maximized the score

T

S

AAAC

AG-C

- We now trace back a path that corresponds to the best alignment

T

S

AAAC

-AGC

AAAC

AG-C

- Sometimes, more than one alignment has the best score

AAAC

A-GC

T

S

Space: O(mn)

Time: O(mn)

- Filling the matrix O(mn)
- Backtrace O(m+n)

- In real-life applications, n and m can be very large
- The space requirements of O(mn) can be too demanding
- If m = n = 1000, we need 1MB space
- If m = n = 10000, we need 100MB space

- We can afford to perform extra computation to save space
- Looping over million operations takes less than seconds on modern workstations

- Can we trade space with time?

A

G

C

1

2

3

0

0

0

-2

-4

-6

1

-2

1

-1

-3

A

2

-4

-1

0

-2

A

3

-6

-3

-2

-1

A

4

-8

-5

-4

-1

C

To compute V[n,m]=d(s[1..n],t[1..m]), we need only O(min(n,m)) space:

- Compute V(i,j), column by column, storing only two columns in memory (or line by line ).

Note however that

- This “trick” fails when we need to reconstruct the optimizing sequence.
- Trace back information requires O(mn) memory bytes.

t

- Construct two alignments
- A= s[1,n/2] vs t[1,j]
- B= s[n/2+1,n] vs t[j+1,m]

s

- Return the concatenated alignment AB

Input: Sequences s[1,n] and t[1,m] to be aligned.

Idea: perform divide and conquer

- Find position (n/2, j) at which the best alignment crosses a midpoint

- We need to find j that maximizes this score. Such j determines a point (n/2,j) through which the best alignment passes.

t

- Thus, we need to compute these two quantities for all values of j

s

- The score of the best alignment that goes through j equals: d(s[1,n/2],t[1,j]) + d(s[n/2+1,n],t[j+1,m])

t

s

Define

- V[i,j] = d(s[1..i],t[1..j]) (As before)
- B[i,j] = d(s[i+1..n],t[j+1..m]) (Symmetrically)
- F[i,j] + B[i,j] = score of best alignment through (i,j)

V[i,j] = d(s[1..i],t[1..j])

As before:

Requires linear space complexity

B[i,j] = d(s[i+1..n],t[j+1..m])

Symmetrically (replacing i with i+1, and j with j+1):

Requires linear space complexity

t

s

- Time to find a mid-point: cnm (c - a constant)
- Size of recursive sub-problems is (n/2,j) and (n/2,m-j-1), hence
T(n,m) = cnm + T(n/2,j) + T(n/2,m-j-1)

Lemma: T(n,m) 2cnm

Proof (by induction):

T(n,m) cnm + 2c(n/2)j + 2c(n/2)(m-j-1) 2cnm.

Thus, time complexity is linear in size of the problem.

At worse, twice the cost of the regular solution.

Consider now a different question:

- Can we find similar substrings of s and t
- Formally, given s[1..n] and t[1..m] find i,j,k, and l such that d(s[i..j],t[k..l]) is maximal

- As before, we use dynamic programming
- We now want to setV[i,j] to record the best alignment of a suffix of s[1..i] and a suffix of t[1..j]
- How should we change the recurrence rule?
- Same as before but with an option to start afresh

- The result is called the Smith-Waterman algorithm

New option:

- We can start a new match instead of extending a previous alignment

Alignment of empty suffixes

T

S

s =TAATA

t =TACTAA

T

S

s =TAATA

t =TACTAA

T

S

s =TAATA

t =TACTAA

T

S

s =TAATA

t =TACTAA

A maximal alignment starts at a maximum entry and follows backwards till a zero entry.

T

S

s =TAATA

t =TACTAA

We have seen two variants of sequence alignment:

- Global alignment
- Local alignment
Other variants in the book and in tutorial time:

- Finding best overlap
- Using an affine cost d(g) = -d –(g-1)e for gaps of length g. The –d is for introducing a gap and –e for continuing the gap. We used d=e=2. We could use smaller e.
These variants are based on the same basic idea of dynamic programming.

Instead of speaking about the score of an alignment, one often talks about an edit distance between two sequences, defined to be the “cost” of the “cheapest” set of edit operations needed to transform one sequence into the other.

- Cheapest operation is “no change”
- Next cheapest operation is “replace”
- The most expensive operation is “add space”.
Our goal is now to minimize the cost of operations, which is equivalent to maximizing the corresponding score (e.g., change scores to edit costs as follows: 1 -1,-1 1,-22).