330 likes | 453 Views
Sequence Alignment. Evolution. Evolution at the DNA level. Deletion. Mutation. …AC GGTG CAGT T ACCA…. SEQUENCE EDITS. …AC ---- CAGT C CACCA…. REARRANGEMENTS. Inversion. Translocation. Duplication. Sequence conservation implies function. Alignment is the key to
E N D
Evolution at the DNA level Deletion Mutation …ACGGTGCAGTTACCA… SEQUENCE EDITS …AC----CAGTCCACCA… REARRANGEMENTS Inversion Translocation Duplication
Sequence conservation implies function • Alignment is the key to • Finding important regions • Determining function • Uncovering the evolutionary forces
Sequence Alignment AGGCTATCACCTGACCTCCAGGCCGATGCCC TAGCTATCACGACCGCGGTCGATTTGCCCGAC -AGGCTATCACCTGACCTCCAGGCCGA--TGCCC--- TAG-CTATCAC--GACCGC--GGTCGATTTGCCCGAC Definition Given two strings x = x1x2...xM, y = y1y2…yN, an alignment is an assignment of gaps to positions 0,…, N in x, and 0,…, N in y, so as to line up each letter in one sequence with either a letter, or a gap in the other sequence
What is a good alignment? AGGCTAGTT, AGCGAAGTTT AGGCTAGTT- 6 matches, 3 mismatches, 1 gap AGCGAAGTTT AGGCTA-GTT- 7 matches, 1 mismatch, 3 gaps AG-CGAAGTTT AGGC-TA-GTT- 7 matches, 0 mismatches, 5 gaps AG-CG-AAGTTT
Scoring Function • Sequence edits: AGGCCTC • Mutations AGGACTC . • Insertions AGGGCCTC • Deletions AGG . CTC Scoring Function: Match: +m Mismatch: -s Gap: -d Score F = (# matches) m - (# mismatches) s – (#gaps) d
How do we compute the best alignment? AGTGCCCTGGAACCCTGACGGTGGGTCACAAAACTTCTGGA Too many possible alignments: >> 2N AGTGACCTGGGAAGACCCTGACCCTGGGTCACAAAACTC
Alignment is additive Observation: The score of aligning x1……xM y1……yN is additive Say that x1…xi xi+1…xM aligns to y1…yj yj+1…yN The two scores add up: F(x[1:M], y[1:N]) = F(x[1:i], y[1:j]) + F(x[i+1:M], y[j+1:N])
Dynamic Programming • There are only a polynomial number of subproblems • Align x1…xi to y1…yj • Original problem is one of the subproblems • Align x1…xM to y1…yN • Each subproblem is easily solved from smaller subproblems • ??? • Then, we can apply Dynamic Programming!!! Let F(i,j) = optimal score of aligning x1……xi y1……yj
Dynamic Programming (cont’d) Notice three possible cases: • xi aligns to yj x1……xi-1 xi y1……yj-1 yj 2. xi aligns to a gap x1……xi-1 xi y1……yj - • yj aligns to a gap x1……xi - y1……yj-1 yj m, if xi = yj F(i,j) = F(i-1, j-1) + -s, if not F(i,j) = F(i-1, j) - d F(i,j) = F(i, j-1) - d
Dynamic Programming (cont’d) How do we know which case is correct? Inductive assumption: F(i, j-1), F(i-1, j), F(i-1, j-1) are optimal Then, F(i-1, j-1) + s(xi, yj) F(i, j) = max F(i-1, j) – d F( i, j-1) – d Where s(xi, yj) = m, if xi = yj; -s, if not
Pairwise Sequence Alignment • As we’ve seen, sequence similarity is an indicator of homology • There are other uses for sequence similarity • Database queries • Comparative genomics • …
Pairwise Sequence Alignment • Example • Which one is better? HEAGAWGHEE PAWHEAE HEAGAWGHE-E HEAGAWGHE-E P-A--W-HEAE --P-AW-HEAE
Scoring • To compare two sequence alignments, calculate a score • PAM or BLOSUM matrices • Matches and mismatches • Gap penalty • Initiating a gap • Gap extension penalty • Extending a gap
C W W -8 17 PAM 250 A R N D C Q E G H I L K M F P S T W Y V B Z A 2 -2 0 0 -2 0 0 1 -1 -1 -2 -1 -1 -3 1 1 1 -6 -3 0 2 1 R -2 6 0 -1 -4 1 -1 -3 2 -2 -3 3 0 -4 0 0 -1 2 -4 -2 1 2 N 0 0 2 2 -4 1 1 0 2 -2 -3 1 -2 -3 0 1 0 -4 -2 -2 4 3 D 0 -1 2 4 -5 2 3 1 1 -2 -4 0 -3 -6 -1 0 0 -7 -4 -2 5 4 C -2 -4 -4 -5 12 -5 -5 -3 -3 -2 -6 -5 -5 -4 -3 0 -2 -8 0 -2 -3 -4 Q 0 1 1 2 -5 4 2 -1 3 -2 -2 1 -1 -5 0 -1 -1 -5 -4 -2 3 5 E 0 -1 1 3 -5 2 4 0 1 -2 -3 0 -2 -5 -1 0 0 -7 -4 -2 4 5 G 1 -3 0 1 -3 -1 0 5 -2 -3 -4 -2 -3 -5 0 1 0 -7 -5 -1 2 1 H -1 2 2 1 -3 3 1 -2 6 -2 -2 0 -2 -2 0 -1 -1 -3 0 -2 3 3 I -1 -2 -2 -2 -2 -2 -2 -3 -2 5 2 -2 2 1 -2 -1 0 -5 -1 4 -1 -1 L -2 -3 -3 -4 -6 -2 -3 -4 -2 2 6 -3 4 2 -3 -3 -2 -2 -1 2 -2 -1 K -1 3 1 0 -5 1 0 -2 0 -2 -3 5 0 -5 -1 0 0 -3 -4 -2 2 2 M -1 0 -2 -3 -5 -1 -2 -3 -2 2 4 0 6 0 -2 -2 -1 -4 -2 2 -1 0 F -3 -4 -3 -6 -4 -5 -5 -5 -2 1 2 -5 0 9 -5 -3 -3 0 7 -1 -3 -4 P 1 0 0 -1 -3 0 -1 0 0 -2 -3 -1 -2 -5 6 1 0 -6 -5 -1 1 1 S 1 0 1 0 0 -1 0 1 -1 -1 -3 0 -2 -3 1 2 1 -2 -3 -1 2 1 T 1 -1 0 0 -2 -1 0 0 -1 0 -2 0 -1 -3 0 1 3 -5 -3 0 2 1 W -6 2 -4 -7 -8 -5 -7 -7 -3 -5 -2 -3 -4 0 -6 -2 -5 17 0 -6 -4 -4 Y -3 -4 -2 -4 0 -4 -4 -5 0 -1 -1 -4 -2 7 -5 -3 -3 0 10 -2 -2 -3 V 0 -2 -2 -2 -2 -2 -2 -1 -2 4 2 -2 2 -1 -1 -1 0 -6 -2 4 0 0 B 2 1 4 5 -3 3 4 2 3 -1 -2 2 -1 -3 1 2 2 -4 -2 0 6 5 Z 1 2 3 4 -4 5 5 1 3 -1 -1 2 0 -4 1 1 1 -4 -3 0 5 6
Example • Gap penalty: -8 • Gap extension: -8 HEAGAWGHE-E --P-AW-HEAE (-8) + (-8) + (-1) + 5 + 15 + (-8) + 10 + 6 + (-8) + 6 = 9 HEAGAWGHE-E Exercise: Calculate for P-A--W-HEAE
Formal Description • Problem:PairSeqAlign • Input: Two sequences x,y Scoring matrix s Gap penalty d Gap extension penalty e • Output: The optimal sequence alignment
How Difficult Is This? • Consider two sequences of length n • There are possible global alignments, and we need to find an optimal one from amongst those!
So what? • So at n = 20, we have over 120 billion possible alignments • We want to be able to align much, much longer sequences • Some proteins have 1000 amino acids • Genes can have several thousand base pairs
Dynamic Programming • General algorithmic development technique • Reuses the results of previous computations • Store intermediate results in a table for reuse • Look up in table for earlier result to build from
Global Alignment • Needleman-Wunsch 1970 • Idea: Build up optimal alignment from optimal alignments of subsequences HEAG --P- -25 Add score from table HEAG- --P-A -33 HEAGA --P-A -20 HEAGA --P— -33 Gap with bottom Top and bottom Gap with top
Global Alignment • Notation • xi – ith letter of string x • yj – jth letter of string y • x1..i – Prefix of x from letters 1 through I • F – matrix of optimal scores • F(i,j) represents optimal score lining up x1..i with y1..j • d – gap penalty • s – scoring matrix MSCS 230: Bioinformatics I - Pairwise Sequence Alignment
Global Alignment • The work is to build up F • Initialize: F(0,0) = 0, F(i,0) = id, F(0,j)=jd • Fill from top left to bottom right using the recursive relation MSCS 230: Bioinformatics I - Pairwise Sequence Alignment
Global Alignment yj aligned to gap Move ahead in both s(xi,yj) d d xi aligned to gap While building the table, keep track of where optimal score came from, reverse arrows MSCS 230: Bioinformatics I - Pairwise Sequence Alignment
The Needleman-Wunsch Matrix x1 ……………………………… xM Every nondecreasing path from (0,0) to (M, N) corresponds to an alignment of the two sequences y1 ……………………………… yN An optimal alignment is composed of optimal subalignments
Performance • Time: O(NM) • Space: O(NM)
Design a Perl Program • Is it doable by Perl? • How can we handle two-dimensional array?
i-> j | V M S M[i, j] = S[j*len+i]