Sequence alignment i lecture 2
This presentation is the property of its rightful owner.
Sponsored Links
1 / 44

Sequence Alignment I Lecture #2 PowerPoint PPT Presentation


  • 84 Views
  • Uploaded on
  • Presentation posted in: General

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

Download Presentation

Sequence Alignment I Lecture #2

An Image/Link below is provided (as is) to download presentation

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 - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript


Sequence alignment i lecture 2

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.

.


Sequence comparison

Sequence Comparison

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


Sequence comparison cont

g1

g2

Sequence Comparison (cont)

  • 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


Sequence alignment

Sequence Alignment

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


Alignments

Alignments

-GCGC-ATGGATTGAGCGA

TGCGCCATTGAT-GACC-A

Three elements:

  • Perfect matches

  • Mismatches

  • Insertions & deletions (indel)


Choosing alignments

Choosing Alignments

There are many possible alignments

For example, compare:

-GCGC-ATGGATTGAGCGA

TGCGCCATTGAT-GACC-A

to

------GCGCATGGATTGAGCGA

TGCGCC----ATTGATGACCA--

Which one is better?


Scoring alignments

Scoring Alignments

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


Simple scoring rule

Simple Scoring Rule

Score each position independently:

  • Match: +1

  • Mismatch: -1

  • Indel -2

    Score of an alignment is sum of position scores


Example

Example

Example:

-GCGC-ATGGATTGAGCGA

TGCGCCATTGAT-GACC-A

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

------GCGCATGGATTGAGCGA

TGCGCC----ATTGATGACCA--

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


More general scores

More General Scores

  • 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

  • Probabilistic interpretation: (e.g.) How likely is one alignment versus another ?


Additive scoring rules

Additive Scoring Rules

  • 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 score

The Optimal Score

  • 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.


Computing optimal score

Computing Optimal Score

  • 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.


Recursive argument

Recursive Argument

  • 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] )


Recursive argument1

Recursive Argument

  • 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] )


Recursive argument2

Recursive Argument

  • 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] )


Recursive argument3

Recursive Argument

Define the notation:

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


Recursive argument4

T

S

Recursive Argument

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

AA

- -

versus

We fill the matrix using the recurrence rule:


Dynamic programming algorithm

T

S

Dynamic Programming Algorithm

We continue to fill the matrix using the recurrence rule


Dynamic programming algorithm1

T

S

+1

-2 -A

A-

-2 (A- versus -A)

Dynamic Programming Algorithm

versus


Dynamic programming algorithm2

T

S

Dynamic Programming Algorithm


Dynamic programming algorithm3

T

S

Conclusion: d(AAAC,AGC) = -1

Dynamic Programming Algorithm


Reconstructing the best alignment

T

S

Reconstructing the Best Alignment

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


Reconstructing the best alignment1

T

S

AAAC

AG-C

Reconstructing the Best Alignment

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


Reconstructing the best alignment2

T

S

AAAC

-AGC

AAAC

AG-C

Reconstructing the Best Alignment

  • Sometimes, more than one alignment has the best score

AAAC

A-GC


Time complexity

T

S

Time Complexity

Space: O(mn)

Time: O(mn)

  • Filling the matrix O(mn)

  • Backtrace O(m+n)


Space complexity

Space Complexity

  • 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?


Why do we need so much space

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

Why Do We Need So Much Space?

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.


Space efficient version outline

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

Space Efficient Version: Outline

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


Finding the 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

Finding the Midpoint

  • 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])


Finding the midpoint algorithm

t

s

Finding the Midpoint (Algorithm)

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)


Computing v i j

Computing V[i,j]

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

As before:

Requires linear space complexity


Computing b i j

Computing B[i,j]

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


Time complexity analysis

t

s

Time Complexity Analysis

  • 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.


Local alignment

Local Alignment

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


Local alignment1

Local Alignment

  • 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


Local alignment2

Local Alignment

New option:

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

Alignment of empty suffixes


Local alignment example

T

S

Local Alignment Example

s =TAATA

t =TACTAA


Local alignment example1

T

S

Local Alignment Example

s =TAATA

t =TACTAA


Local alignment example2

T

S

Local Alignment Example

s =TAATA

t =TACTAA


Local alignment example3

T

S

Local Alignment Example

s =TAATA

t =TACTAA

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


Local alignment example4

T

S

Local Alignment Example

s =TAATA

t =TACTAA


Variants of sequence alignment

Variants of Sequence Alignment

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.


Remark edit distance

Remark: Edit Distance

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,-22).


  • Login