CS CS 466 466 In Introduct ctio ion t to B Bio ioin - - PowerPoint PPT Presentation

CS CS 466 466 In Introduct ctio ion t to B Bio ioin informatics ics

Lecture 6

Mohammed El-Kebir February 6, 2020

Course Announcements

Instructor:

• Mohammed El-Kebir (melkebir)
• Office hours: Wednesdays, 3:15-4:15pm

TA:

• Aswhin Ramesh (aramesh7)
• Office hours: Fridays, 11:00-11:59am in SC 3405

2

Homework 1: Due on Sept. 18 (11:59pm) Midterm: 10/4, 11-1pm @Transportation Building 103 (conflict: 10/7, 7-9pm @Siebel 1302 -- to sign up email me)

Global, Fitting and Local Alignment

3

Local Alignment problem: Given strings 𝐰 ∈ Σ$ and 𝐱 ∈ Σ& and scoring function 𝜀, find a substring

• f 𝐰 and a substring of 𝐱 whose alignment has

maximum global alignment score 𝑡∗ among all global alignments of all substrings of 𝐰 and 𝐱 [Smith-Waterman algorithm] Fitting Alignment problem: Given strings 𝐰 ∈ Σ$ and 𝐱 ∈ Σ& and scoring function 𝜀, find an alignment of 𝐰 and a substring of 𝐱 with maximum global alignment score 𝑡∗ among all global alignments of 𝐰 and all substrings of 𝐱 Global Alignment problem: Given strings 𝐰 ∈ Σ$ and 𝐱 ∈ Σ& and scoring function 𝜀, find alignment

• f 𝐰 and 𝐱 with maximum score.

[Needleman-Wunsch algorithm]

A G G T A C G G C

𝐰\𝐱

Question: How to assess resulting algorithms?

Time Complexity

4

1 2 3 n = 4 O O O O O 1 O O O O O 2 O O O O O 3 O O O O O m = 4 O O O O O

W A T C G A T G T V Alignment is a path from source (0, 0) to target (𝑛, 𝑜) in edit graph Edit graph is a weighed, directed grid graph 𝐻 = (𝑊, 𝐹) with source vertex (0, 0) and target vertex (𝑛, 𝑜). Each edge ((𝑗, 𝑘), (𝑙, 𝑚)) has weight depending on direction. Running time is 𝑃(𝑛𝑜) [quadratic time]

Time Complexity

5

1 2 3 n = 4 O O O O O 1 O O O O O 2 O O O O O 3 O O O O O m = 4 O O O O O

W A T C G A T G T V Alignment is a path from source (0, 0) to target (𝑛, 𝑜) in edit graph Edit graph is a weighed, directed grid graph 𝐻 = (𝑊, 𝐹) with source vertex (0, 0) and target vertex (𝑛, 𝑜). Each edge ((𝑗, 𝑘), (𝑙, 𝑚)) has weight depending on direction. Running time is 𝑃(𝑛𝑜) [quadratic time] Question: Compute alignment faster than 𝑃(𝑛𝑜) time? [subquadratic time]

Space Complexity

6

1 2 3 n = 4 1 2 3 m = 4

W A T C G A T G T V Thus, space complexity is 𝑃(𝑛𝑜) [quadratic space] Size of DP table is 𝑛 + 1 × (𝑜 + 1) Example: To align a short read (𝑛 = 100) to human genome (𝑜 = 3 = 10>), we need 300 GB memory.

Space Complexity

7

1 2 3 n = 4 1 2 3 m = 4

W A T C G A T G T V Thus, space complexity is 𝑃(𝑛𝑜) [quadratic space] Size of DP table is 𝑛 + 1 × (𝑜 + 1) Example: To align a short read (𝑛 = 100) to human genome (𝑜 = 3 = 10>), we need 300 GB memory. Question: How long is an alignment?

Space Complexity

8

1 2 3 n = 4 1 2 3 m = 4

W A T C G A T G T V Thus, space complexity is 𝑃(𝑛𝑜) [quadratic space] Size of DP table is 𝑛 + 1 × (𝑜 + 1) Example: To align a short read (𝑛 = 100) to human genome (𝑜 = 3 = 10>), we need 300 GB memory. Question: How long is an alignment? Question: Compute alignment in 𝑃(𝑛) space? [linear space]

Outline

• 1. Recap of global, fitting, local and gapped alignment
• 2. Space-efficient alignment
• 3. Subquadratic time alignment

Reading:

• Jones and Pevzner. Chapters 7.1-7.4
• Lecture notes

9

𝑛 𝑜

Space Efficient Alignment

10

Computing 𝑡[𝑗, 𝑘] requires access to: 𝑡 𝑗 − 1, 𝑘 , 𝑡[𝑗, 𝑘 − 1] and 𝑡[𝑗 − 1, 𝑘 − 1] 𝑗 𝑘

s[i, j] = max          0, if i = 0 and j = 0, s[i − 1, j] + δ(vi, −), if i > 0, s[i, j − 1] + δ(−, wj), if j > 0, s[i − 1, j − 1] + δ(vi, wj), if i > 0 and j > 0.

𝑛 𝑜

Space Efficient Alignment

11

Computing 𝑡[𝑗, 𝑘] requires access to: 𝑡 𝑗 − 1, 𝑘 , 𝑡[𝑗, 𝑘 − 1] and 𝑡[𝑗 − 1, 𝑘 − 1] 𝑗 𝑘 Thus it suffices to store only two columns to compute optimal alignment score 𝑡 𝑛, 𝑜 , i.e., 2 𝑛 + 1 = 𝑃(𝑛) space.

s[i, j] = max          0, if i = 0 and j = 0, s[i − 1, j] + δ(vi, −), if i > 0, s[i, j − 1] + δ(−, wj), if j > 0, s[i − 1, j − 1] + δ(vi, wj), if i > 0 and j > 0.

𝑛 𝑜

Space Efficient Alignment

12

Question: What if we want alignment itself? Computing 𝑡[𝑗, 𝑘] requires access to: 𝑡 𝑗 − 1, 𝑘 , 𝑡[𝑗, 𝑘 − 1] and 𝑡[𝑗 − 1, 𝑘 − 1] 𝑗 𝑘 Thus it suffices to store only two columns to compute optimal alignment score 𝑡 𝑛, 𝑜 , i.e., 2 𝑛 + 1 = 𝑃(𝑛) space.

s[i, j] = max          0, if i = 0 and j = 0, s[i − 1, j] + δ(vi, −), if i > 0, s[i, j − 1] + δ(−, wj), if j > 0, s[i − 1, j − 1] + δ(vi, wj), if i > 0 and j > 0.

Space Efficient Alignment – First Attempt

• What if also want optimal alignment?
• Easy: keep best pointers as fill in table.
• No! Do not know which path to keep until computing

recurrence at each step.

w v w v

Space Efficient Alignment – First Attempt

• What if also want optimal alignment?
• Easy: keep best pointers as fill in table.
• No! Do not know which path to keep until computing

recurrence at each step.

w v w v

Space Efficient Alignment – First Attempt

• What if also want optimal alignment?
• Easy: keep best pointers as fill in table.
• No! Do not know which path to keep until computing

recurrence at each step.

w v w v

Best score for column might not be part of best alignment!

Space Efficient Alignment – Second Attempt

16

𝑜/2 𝑛 𝑗∗ 𝑗

Maximum weight path from (0,0) to (𝑛, 𝑜) passes through (𝑗∗, 𝑜/2) Question: What is 𝑗∗? Alignment is a path from source (0, 0) to target (𝑛, 𝑜) in edit graph

17

Hirschberg(𝑗, 𝑘, 𝑗E, 𝑘′) 1. if 𝑘E − 𝑘 > 1 2. 𝑗∗ ß arg max

MNMOONMO wt(𝑗′′)

3. Report (𝑗∗, 𝑘 + ROSR

T )

4. Hirschberg(𝑗, 𝑘, 𝑗∗, 𝑘 + ROSR

T )

5. Hirschberg(𝑗∗, 𝑘 + ROSR

T , 𝑗E, 𝑘′)

18

Time: area + area/2 + area/4 + … = area (1 + ½ + ¼ + ⅛ + …) ≤ 2 × area = O(mn) Space: O(m) Hirschberg(𝑗, 𝑘, 𝑗E, 𝑘′) 1. if 𝑘E − 𝑘 > 1 2. 𝑗∗ ß arg max

MNMOONMO wt(𝑗′′)

3. Report (𝑗∗, 𝑘 + ROSR

T )

4. Hirschberg(𝑗, 𝑘, 𝑗∗, 𝑘 + ROSR

T )

5. Hirschberg(𝑗∗, 𝑘 + ROSR

T , 𝑗E, 𝑘′)

19

Time: area + area/2 + area/4 + … = area (1 + ½ + ¼ + ⅛ + …) ≤ 2 × area = O(mn) Space: O(m) Question: How to reconstruct alignment from reported vertices? Hirschberg(𝑗, 𝑘, 𝑗E, 𝑘′) 1. if 𝑘E − 𝑘 > 1 2. 𝑗∗ ß arg max

MNMOONMO wt(𝑗′′)

3. Report (𝑗∗, 𝑘 + ROSR

T )

4. Hirschberg(𝑗, 𝑘, 𝑗∗, 𝑘 + ROSR

T )

5. Hirschberg(𝑗∗, 𝑘 + ROSR

T , 𝑗E, 𝑘′)

Hirschberg Algorithm: Reversing Edges Necessary?

20

𝑗∗ = arg max{

VNMN$

preYix 𝑗 + sufYix 𝑗 } Max weight path from (0,0) to (𝑛, 𝑜) through (𝑗∗, 𝑜/2) Compute preYix 𝑗 0 ≤ 𝑗 ≤ 𝑛} in O(𝑛𝑘) time and O(𝑛) space, by starting from (0,0) to 𝑛, 𝑘 keeping only two columns in memory. [single-source multiple destinations]

𝑘 𝑛 𝑗∗ 𝑗

Hirschberg Algorithm: Reversing Edges Necessary?

21

𝑗∗ = arg max{

VNMN$

preYix 𝑗 + sufYix 𝑗 } Max weight path from (0,0) to (𝑛, 𝑜) through (𝑗∗, 𝑜/2) Want: Compute sufYix 𝑗 0 ≤ 𝑗 ≤ 𝑛} in O(𝑛𝑘) time and O(𝑛) space Doing a longest path from each 𝑗, 𝑘 to 𝑛, 𝑜 (for all 0 ≤ 𝑗 ≤ 𝑛) will not achieve desired running time! Reversing edges enables single-source multiple destination computation in desired time and space bound!

𝑘 𝑛 𝑗∗ 𝑗

Compute preYix 𝑗 0 ≤ 𝑗 ≤ 𝑛} in O(𝑛𝑘) time and O(𝑛) space, by starting from (0,0) to 𝑛, 𝑘 keeping only two columns in memory. [single-source multiple destinations]

Hirschberg Algorithm: Reconstructing Alignment

22

Problem: Given reported vertices and scores { 𝑗V, 0, 𝑡V , … , 𝑗&, 𝑜, 𝑡& }, find intermediary vertices. Hirschberg(𝑗, 𝑘, 𝑗E, 𝑘′) 1. if 𝑘E − 𝑘 > 1 2. 𝑗∗ ß arg max

VNMN$

wt(𝑗) 3. Report (𝑗∗, 𝑘 + ROSR

T )

4. Hirschberg(𝑗, 𝑘, 𝑗∗, 𝑘 + ROSR

T )

5. Hirschberg(𝑗∗, 𝑘 + ROSR

T , 𝑗E, 𝑘′)

1 2 3 4 5

• 1
• 2
• 3
• 4
• 5

1

• 1

1

• 1
• 2
• 3

2

• 2

2 1

• 1

3

• 3
• 1

1 1 2 1 4

• 4
• 2

1 1 5

• 5
• 3
• 1

1 2

W A T C G A T G T V A T

• G

T C A T C G

• C

C C

Transposing matrix does not help, because gaps could occur in both input sequences

Linear Space Alignment – The Hirschberg Algorithm

23

Outline

• 1. Recap of global, fitting, local and gapped alignment
• 2. Space-efficient alignment
• 3. Subquadratic time alignment

Reading:

• Jones and Pevzner. Chapters 7.1-7.4
• Lecture notes

24

Banded Alignment

25

x1 x2 x3 . . . . . . xM y1 y2 y3 ... ... ... ... yN

Constrain traceback to band of DP matrix (penalize big gaps)

Figure source: http://jinome.stanford.edu/stat366/pdfs/stat366_win0607_lecture04.pdf

Constraint path to band of width 𝑙 around diagonal Running time: O(𝑜𝑙) Gives a good approximation of highly identical sequences

Banded Alignment

26

x1 x2 x3 . . . . . . xM y1 y2 y3 ... ... ... ... yN

Constrain traceback to band of DP matrix (penalize big gaps)

Figure source: http://jinome.stanford.edu/stat366/pdfs/stat366_win0607_lecture04.pdf

Constraint path to band of width 𝑙 around diagonal Running time: O(𝑜𝑙) Gives a good approximation of highly identical sequences Question: How to change recurrence to accomplish this?

Block Alignment

27

Divide input sequences into blocks of length 𝑢

v1, …, vt vt+1, …, v2t … vm-t+1, …, vm w1, …, wt wt+1, …, w2t … wn-t+1, …, wn

Block Alignment

28

Divide input sequences into blocks of length 𝑢

v1, …, vt vt+1, …, v2t … vm-t+1, …, vm w1, …, wt wt+1, …, w2t … wn-t+1, …, wn

Require that paths in edit graph pass through corners of blocks

Block Alignment

29

Divide input sequences into blocks of length 𝑢

v1, …, vt vt+1, …, v2t … vm-t+1, …, vm w1, …, wt wt+1, …, w2t … wn-t+1, …, wn

Require that paths in edit graph pass through corners of blocks

s[i, j] = max          0, if i = 0 and j = 0, s[i − 1, j] − σ, if i > 0, s[i, j − 1] − σ, if j > 0, s[i − 1, j − 1] + β(i, j), if i > 0 and j > 0.

0 ≤ 𝑗, 𝑘 ≤ 𝑢 and 𝛾(𝑗, 𝑘) is max score alignment between block 𝑗 of 𝐰 and block 𝑘 of 𝐱

Block Alignment – First Attempt: Pre-compute 𝛾 𝑗, 𝑘

30

s[i, j] = max          0, if i = 0 and j = 0, s[i − 1, j] − σ, if i > 0, s[i, j − 1] − σ, if j > 0, s[i − 1, j − 1] + β(i, j), if i > 0 and j > 0.

0 ≤ 𝑗, 𝑘 ≤ 𝑜/𝑢 and 𝛾(𝑗, 𝑘) is max score alignment between block 𝑗 of 𝐰 and block 𝑘 of 𝐱

t 2t nt t 2t nt

Block Alignment – First Attempt: Pre-compute 𝛾 𝑗, 𝑘

31

s[i, j] = max          0, if i = 0 and j = 0, s[i − 1, j] − σ, if i > 0, s[i, j − 1] − σ, if j > 0, s[i − 1, j − 1] + β(i, j), if i > 0 and j > 0.

0 ≤ 𝑗, 𝑘 ≤ 𝑜/𝑢 and 𝛾(𝑗, 𝑘) is max score alignment between block 𝑗 of 𝐰 and block 𝑘 of 𝐱 Question: How much time to compute all 𝛾(𝑗, 𝑘)?

t 2t nt t 2t nt

Block Alignment – First Attempt: Pre-compute 𝛾 𝑗, 𝑘

32

s[i, j] = max          0, if i = 0 and j = 0, s[i − 1, j] − σ, if i > 0, s[i, j − 1] − σ, if j > 0, s[i − 1, j − 1] + β(i, j), if i > 0 and j > 0.

0 ≤ 𝑗, 𝑘 ≤ 𝑜/𝑢 and 𝛾(𝑗, 𝑘) is max score alignment between block 𝑗 of 𝐰 and block 𝑘 of 𝐱 Question: How much time to compute all 𝛾(𝑗, 𝑘)?

t 2t nt t 2t nt

Computing 𝛾 𝑗, 𝑘 takes 𝑃 𝑢T time There are 𝑜/𝑢 × 𝑜/𝑢 values 𝛾 𝑗, 𝑘 Total: 𝑃

& d × & d × 𝑢T = 𝑃 𝑜T time

Block Alignment – Four Russians Technique

33

t 2t nt t 2t nt

Pre-compute and store all βij Pre-compute and store all max weight alignments 𝑇[𝐰′, 𝐱′] of all pairs (𝐰E, 𝐱E) of length t strings

Algorithm:

• 1. Precompute 𝑇[𝐰′, 𝐱′] where

𝐰E, 𝐱E ∈ Σd

• 2. Compute block alignment

between 𝐰 and 𝐱 using 𝑇

s[i, j] = max          0, if i = 0 and j = 0, s[i − 1, j] − σ, if i > 0, s[i, j − 1] − σ, if j > 0, s[i − 1, j − 1] + S[v(i), w(j)], if i > 0 and j > 0.

Block Alignment – Four Russians Technique

34

t 2t nt t 2t nt

Pre-compute and store all βij Pre-compute and store all max weight alignments 𝑇[𝐰′, 𝐱′] of all pairs (𝐰E, 𝐱E) of length t strings

Algorithm:

• 1. Precompute 𝑇[𝐰′, 𝐱′] where

𝐰E, 𝐱E ∈ Σd

• 2. Compute block alignment

between 𝐰 and 𝐱 using 𝑇

s[i, j] = max          0, if i = 0 and j = 0, s[i − 1, j] − σ, if i > 0, s[i, j − 1] − σ, if j > 0, s[i − 1, j − 1] + S[v(i), w(j)], if i > 0 and j > 0.

Question: How to choose 𝑢 for DNA?

Fastest Subquadratic Alignment* Algorithm

35

*for edit distance

Edit distance in O(𝑜T/ log 𝑜) time Barely subquadratic! Want: O(𝑜TSh) time where 𝜁 > 0

Fastest Subquadratic Alignment* Algorithm

36

*for edit distance

Edit distance in O(𝑜T/ log 𝑜) time Barely subquadratic! Want: O(𝑜TSh) time where 𝜁 > 0 Question: Is 𝑜TSh in O(𝑜T/ log 𝑜) for any 𝜁 > 0?

Hardness Result for Edit Distance [STOC 2015]

37

38

Take Home Messages

• 1. Global alignment in O(𝑛𝑜) time and O(𝑛) space
• Hirschberg algorithm
• 2. Block alignment can be done in subquadratic time
• Four Russians Technique: O(𝑜T/ log 𝑜) time
• 3. Global alignment cannot be done in O(𝑜TSh) time under SETH

Reading:

• Jones and Pevzner. Chapters 7.1-7.4
• Lecture notes

39