Week2n3_AlignmentDynamicProgramming

Recap of last week (con’d): I/O of a mutation calling tool • Input: 1) a reference genome, and 2) the sequences of an individual’s DNA, which is not exactly the same as the reference genome. • Output: the “equivalent” regions on the reference for each read to maximize their similarity. Reference Individual’s Genome Fragments of individual’s genome

Recap of last week (con’d): sequence length in alignment • Suppose V is the matching score, m and p are the sequence lengths for X and Y . • The final alignment score S is V – m – p. ACGGGTTTT ACGTGTTTT 8 ATGGGAAAATTTTAGCTG ACGTGAGACTCATCGGAG 9 v.s. V = S = 8 – 9 – 9 = -10 9 – 18 – 18 = -27 9 nucleotides 18 nucleotides 9 nucleotides 18 nucleotides >

Materials for the next two weeks • In the last week, we focused on the evaluation of the alignments , suppose we already have the alignments. But we didn’t discuss how to come up with the alignment originally. • Today we will see how computers can be taught to come up with the alignments automatically . • To make this happen, we need to compose the algorithms (a process or set of rules to be followed in calculations or other problem-solving operations). • We will focus on one algorithm for alignment, called dynamic programming .

Why shall we bother to come up with an algorithm? • A naïve method is to score all alignments, and select the one with the highest score. • But there will be exponential number of alignments to score (depending on ࠵? and ࠵? ), and it is unrealistic to calculate all of them when ࠵? and ࠵? are large. • The dynamic programming algorithm saves us tremendous amount of computing time (and computer memory) for finding the highest- scoring alignment.

What does it mean? Let’s see an example. • Suppose we know • the optimal alignment between ACGACT and ACT is and the optimal score is F(ACGACT, ACT) • the optimal alignment between ACGACT T and ACT is and the optimal score is F(ACGACT T , ACT) • the optimal alignment between ACGACT and ACT T is and the optimal score is F(ACGACT, ACT T ) ACGACT TTG ACT TTCG Optimal score of ACGACT T and ACT T = ? ACGACTT TG ACT TTCG ACGACT TTG ACTT TCG F(ACGACT, ACT) F(ACGACT T , ACT) F(ACGACT, ACT T )

Now let’s generalize it. ACGACT A - - -CT ACGACTT A - - -CT- ACGACT T A - - -CT T ACGACTT - A - - - CT- T ACGACT- A - - -CTT ACGACT- T A - - -CTT - + ࠵?(࠵? ! , ࠵? " ) - d - d 1 2 3 Notation: • ࠵? ! : subsequence btw. the 1 st a nd ࠵? th nucleotide on ࠵? ; • ࠵? " : subsequence btw. the 1 st a nd ࠵? th nucleotide on ࠵? ; Suppose ࠵? is the index of the nucleotide of the first sequence ࠵? , and ࠵? is the index of the nucleotide of the second sequence ࠵? . The optimal alignment score of ACGACTT and ACTT is ࠵?(࠵? ࠵? , ࠵? ࠵? ) , ࠵? = 7 , ࠵? = 4 . ࠵?(࠵? !#$ , ࠵? "#$ ) ࠵? $ = ࠵?(࠵? !#$ , ࠵? "#$ ) ࠵?(࠵? ! , ࠵? "#$ ) ࠵? % = ࠵?(࠵? ! , ࠵? "#$ ) ࠵?(࠵? !#$ , ࠵? " ) ࠵? & = ࠵?(࠵? !#$ , ࠵? " ) ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? !#$ , ࠵? "#$ + ࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵? • ࠵? ! : the ࠵? th nucleotide on ࠵?; • ࠵? " : the ࠵? th nucleotide on ࠵? • ࠵?(࠵? ! , ࠵? " ) : score of aligning ࠵? ! and ࠵? "

Entries necessary to calculate ࠵? ! . ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ࠵? $ ࠵? % ࠵? & ? A C G A C T T T G A C T T T C G ࠵? ࠵?

First of all, what is ࠵?(࠵? ! , ࠵? ! ) ? ? ࠵? $ ࠵? % ࠵? & ? A C G A C T T T G A C T T T C G ࠵? ࠵?

We need an extra row and column ! ? ࠵? $ ࠵? % ࠵? & ? A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵?

How to define/initialize ࠵? ࠵? " , ࠵? # , ࠵? = 1 … ࠵? and ࠵? ࠵? $ , ࠵? " , ࠵? = 1 … ࠵? Depends on which problem we are trying to solve. Global, local or overlap alignment?

Initialization in global alignment

Recap of global alignment • End-to-end matching ACGACTTTG ACTTTCG have a prior knowledge that these two sequences align with each other from beginning to the end Application: • Encourages end-to-end mapping. • For any indels occurring at the beginning or the end of the sequence, penalize it. • Penalty can be given in a linear manner.

Global alignment: 0 ࠵? $ ࠵? % ࠵? & ? A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d

Global alignment: 0 -d ࠵? $ ࠵? % ࠵? & ? A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d

What does it mean by saying ࠵? ࠵? ! , ࠵? " = −࠵? ? It means that if we treat the first nucleotide of ࠵? as being mapped to nowhere on ࠵? (thus a gap), the penalty is −࠵? .

Global alignment: 0 -d ࠵? $ ࠵? % ࠵? & ? A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d linear gap penalty

Global alignment: 0 -d -2d -3d -4d -5d -6d -7d -8d -9d -d -2d -3d ࠵? $ ࠵? % -4d ࠵? & ? -5d -6d -7d A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d linear gap penalty

Suppose we use the following scoring scheme • For indels, ࠵? = 1 • For matches/mismatches: A C G T A C G T 1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 -1 -1 -1 -1 1

Now we will fill in the whole matrix according to 1) initialization of ࠵? ࠵? C , ࠵? D = −࠵?࠵?, ࠵? ࠵? E , ࠵? C = −࠵?࠵? 2) ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? !#$ , ࠵? "#$ + ࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵? ࠵? = 1 … ࠵?, ࠵? = 1 … ࠵?

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 ? -2 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d ࠵? ࠵? $ , ࠵? $ = max @ ࠵? ࠵? ’ , ࠵? ’ + ࠵? ࠵? $ , ࠵? $ = 0 + 1 = 1 ࠵? ࠵? $ , ࠵? ’ − ࠵? = −1 − 1 = −2 ࠵? ࠵? ’ , ࠵? $ − ࠵? = −1 − 1 = −2

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 -2 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d ࠵? ࠵? $ , ࠵? $ = max @ ࠵? ࠵? ’ , ࠵? ’ + ࠵? ࠵? $ , ࠵? $ = 0 + 1 = 1 ࠵? ࠵? $ , ࠵? ’ − ࠵? = −1 − 1 = −2 ࠵? ࠵? ’ , ࠵? $ − ࠵? = −1 − 1 = −2

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 ? -2 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d ࠵? ࠵? % , ࠵? $ = max @ ࠵? ࠵? $ , ࠵? ’ + ࠵? ࠵? % , ࠵? $ = −1 − 1 = −2 ࠵? ࠵? % , ࠵? ’ − ࠵? = −2 − 1 = −3 ࠵? ࠵? $ , ࠵? $ − ࠵? = 1 − 1 = 0

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -2 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d ࠵? ࠵? % , ࠵? $ = max @ ࠵? ࠵? $ , ࠵? ’ + ࠵? ࠵? % , ࠵? $ = −1 − 1 = −2 ࠵? ࠵? % , ࠵? ’ − ࠵? = −2 − 1 = −3 ࠵? ࠵? $ , ࠵? $ − ࠵? = 1 − 1 = 0

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 ? -2 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d ࠵? ࠵? & , ࠵? $ = max @ ࠵? ࠵? % , ࠵? ’ + ࠵? ࠵? & , ࠵? $ = −2 − 1 = −3 ࠵? ࠵? & , ࠵? ’ − ࠵? = −3 − 1 = −4 ࠵? ࠵? % , ࠵? $ − ࠵? = 0 − 1 = −1

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d ࠵? ࠵? & , ࠵? $ = max @ ࠵? ࠵? % , ࠵? ’ + ࠵? ࠵? & , ࠵? $ = −2 − 1 = −3 ࠵? ࠵? & , ࠵? ’ − ࠵? = −3 − 1 = −4 ࠵? ࠵? % , ࠵? $ − ࠵? = 0 − 1 = −1

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 ? -2 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d ࠵? ࠵? ( , ࠵? $ = max @ ࠵? ࠵? & , ࠵? ’ + ࠵? ࠵? ( , ࠵? $ = −3 + 1 = −2 ࠵? ࠵? ( , ࠵? ’ − ࠵? = −4 − 1 = −5 ࠵? ࠵? & , ࠵? $ − ࠵? = −1 − 1 = −2

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -2 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d ࠵? ࠵? ( , ࠵? $ = max @ ࠵? ࠵? & , ࠵? ’ + ࠵? ࠵? ( , ࠵? $ = −3 + 1 = −2 ࠵? ࠵? ( , ࠵? ’ − ࠵? = −4 − 1 = −5 ࠵? ࠵? & , ࠵? $ − ࠵? = −1 − 1 = −2 This will result in two alignments (paths) with the same alignment score.

We will then do the same thing row by row, or column by column, until we fill in the whole matrix.

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 -5 -3 -1 -2 -1 -1 -1 0 2 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d

Now that we have filled the whole matrix, how can we figure out the alignment? We do this by back tracing . In global alignment, we start from the bottom right entry and follow the arrows, but in a backward manner to find the optimal path.

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 -5 -3 -1 -2 -1 -1 -1 0 2 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d Global alignment’s back tracing always start from ࠵?(࠵? # , ࠵? $ ) !

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 -5 -3 -1 -2 -1 -1 -1 0 2 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d We follow the arrows in a backward manner to trace the alignment path.

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 -5 -3 -1 -2 -1 -1 -1 0 2 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d We follow the arrows in a backward manner to trace the alignment path.

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 -5 -3 -1 -2 -1 -1 -1 0 2 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d We follow the arrows in a backward manner to trace the alignment path.

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 -5 -3 -1 -2 -1 -1 -1 0 2 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d We follow the arrows in a backward manner to trace the alignment path.

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 -5 -3 -1 -2 -1 -1 -1 0 2 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d We follow the arrows in a backward manner to trace the alignment path.

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 -5 -3 -1 -2 -1 -1 -1 0 2 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d For multiple paths, we label all of them. They will be equally good alignments.

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 -5 -3 -1 -2 -1 -1 -1 0 2 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d For multiple paths, we label all of them. They will be equally good alignments.

Global alignment: 0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -1 1 0 -1 -2 -3 -4 -5 -6 -7 -2 0 2 1 0 -1 -2 -3 -4 -5 -3 -1 1 1 0 -1 0 -1 -2 -3 -4 -2 0 0 0 -1 0 1 0 -1 -5 -3 -1 -1 -1 -1 0 1 2 1 -6 -4 -2 -2 -2 0 -1 0 1 1 -7 -5 -3 -1 -2 -1 -1 -1 0 2 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In global alignment, we shall back trace it to the upper top entry: ࠵?(࠵? % , ࠵? % ) Notice in this slide, I added the arrows on the initialized entries as well, which has only one choice.

Given all the back tracing, how can we write the alignment between ࠵? and ࠵? ? We do this by putting the matching/mismatching nucleotides corresponding to each other, whereas for the deleted nucleotide, we use symbol “-” to show that the nucleotide is absent at a certain position.

Global alignment: First, we write the sequences on the upper and bottom rows. ACGACTTTG ACTTTCG

Global alignment: First, we write the sequences on the upper and bottom rows. ACGACTTTG ACTTTCG Second, we start from the lower bottom and identify matches/mismatches and indels between the two sequences. a match diagonal: match; horizontal: - on Y; vertical: - on X

Global alignment: First, we write the sequences on the upper and bottom rows. ACGACTTTG ACTTTCG Second, we start from the lower bottom and identify matches/mismatches and indels between the two sequences. a match diagonal: match; horizontal: - on Y; vertical: - on X

Global alignment: First, we write the sequences on the upper and bottom rows. ACGACTTT-G ACTTTCG Second, we start from the lower bottom and identify matches/mismatches and indels between the two sequences. a gap For an indel, add a “-” at the deleted position on the corresponding sequence. diagonal: match; horizontal: - on Y; vertical: - on X

Global alignment: First, we write the sequences on the upper and bottom rows. Second, we start from the lower bottom and identify matches/mismatches and indels between the two sequences. For an indel, add a “-” at the deleted position on the corresponding sequence. We then continue until all the nucleotides find their matches/mismatches or indels. ACGACTTT-G - - - ACTTTCG Optimal alignment #1 diagonal: match; horizontal: - on Y; vertical: - on X

Global alignment: What about this path? Can you write down the alignment? Optimal alignment #2 diagonal: match; horizontal: - on Y; vertical: - on X

Global alignment: ACGACTTT-G A- - -CTTTCG Here is alignment #2. diagonal: match; horizontal: - on Y; vertical: - on X

Global alignment: Can you practice on the 3 rd alignment? Optimal alignment #3 diagonal: match; horizontal: - on Y; vertical: - on X

Global alignment: Here is alignment #3 . ACGACTTT-G AC- - -TTTCG diagonal: match; horizontal: - on Y; vertical: - on X

In total there are three equally optimal global alignment! All have score 2. ACGACTTT-G - - - ACTTTCG ACGACTTT-G AC- - -TTTCG ACGACTTT-G A- - -CTTTCG

The dynamic programming algorithm for global alignment is called Needleman-Wunsch algorithm.

Summary of global alignment (or Needleman- Wunsch algorithm) • Initialization: • Penalize the gaps in the beginning of the two sequences. • In our case, we used linear gap penalty: ࠵? ࠵? ’ , ࠵? " = −࠵?࠵?, ࠵? ࠵? ! , ࠵? ’ = −࠵?࠵? • Recursion or recursive function: • Back tracing: • Start from the lower right entry; • End at the top left entry. ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? !#$ , ࠵? "#$ + ࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Problem definition for local alignment • To align two sequences ࠵? !,…,’ and ࠵? !,…,( finding the best scoring alignment between ࠵? and ࠵?, in which only the subsequences of nucleotides shall be included. • Part of the reads is aligned ACGACTTTG ACTTTCG local alignment

Dynamic programming for pairwise local alignment • Problem definition • Input: sequences ࠵? (of length ࠵? ) and ࠵? (of length ࠵? ) • Output: the alignment between ࠵? and ࠵? involving subsequences of nucleotides for both ࠵? and ࠵? whose alignment score is the highest

How about the local alignment between ࠵? and ࠵? ? A few differences from global alignment: ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? !#$ , ࠵? "#$ + ࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵? ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? ࠵? !#$ , ࠵? "#$ + ࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵? • No penalties for indels in the beginning of the two sequences, i.e., 0’s for all initial entries ࠵? ) , ࠵? = 1, … , ࠵?, ࠵? * , ࠵? = 1, … , ࠵? • The back tracing can start from anywhere, whichever entry has the highest value in the matrix. Back tracing stops at the first zero. • For ࠵?(࠵? ) , ࠵? * ) , instead of: We use:

Local alignment: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In local alignment, we use 0’s for all initializations. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? ࠵? !#$ , ࠵? "#$ +࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Local alignment: 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In local alignment, we use 0’s for all initializations. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? ࠵? !#$ , ࠵? "#$ +࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Local alignment: 0 0 0 0 0 0 0 0 0 0 0 1 ? 0 0 0 0 0 0 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In local alignment, we use 0’s for all initializations. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? ࠵? !#$ , ࠵? "#$ +࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Local alignment: 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In local alignment, we use 0’s for all initializations. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? ࠵? !#$ , ࠵? "#$ +࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Local alignment: 0 0 0 0 0 0 0 0 0 0 0 1 0 ? 0 0 0 0 0 0 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In local alignment, we use 0’s for all initializations. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? ࠵? !#$ , ࠵? "#$ +࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Local alignment: 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In local alignment, we use 0’s for all initializations. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? ࠵? !#$ , ࠵? "#$ +࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵? coming from nowhere, thus no arrows

Now we try to fill in the whole matrix for local alignment.

Local alignment: 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 2 1 0 2 1 0 0 0 0 0 1 1 0 1 3 2 1 0 0 0 0 0 0 0 2 4 3 2 0 0 0 0 0 0 1 3 5 4 0 0 1 0 0 1 0 2 4 4 0 0 0 2 1 0 0 1 3 5 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d Back tracing starts from the entry (entries) with the highest value. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? ࠵? !#$ , ࠵? "#$ +࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Local alignment: ACGACTTT-G ACTTTCG optimal local alignment #1 In back tracing of local alignment, we stop at the first zero. In the corresponding alignment, we don’t penalize the indel in the beginning of the two sequences. ACGACTTTG ACTTTCG optimal local alignment #2

The dynamic programming algorithm for local alignment is called Smith-Waterman algorithm.

Summary of local alignment (or Smith- Waterman algorithm) • Initialization: • No penalty on the gaps in the beginning of the two sequences. • Thus ࠵? ࠵? ’ , ࠵? " = 0, ࠵? ࠵? ! , ࠵? ’ = 0 • Recursion or recursive function: • Back tracing: • Start from the entry/entries with the highest score; • End at the first entry with score zero. ࠵? ࠵? ! , ࠵? " = max 0 ࠵? ࠵? !#$ , ࠵? "#$ + ࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Now, let’s move on to overlap alignment.

Problem definition for overlap alignment • To align two sequences ࠵? !,…,’ and ࠵? !,…,( finding the best scoring alignment between ࠵? and ࠵?, in which either one end for both ࠵? and ࠵? in different directions, or both ends for either ࠵? or ࠵? are involved. • One end of a read aligns to the overhang of another • Part of one read aligns to the other entire read Case #1 Case #2 Case #3 Case #4

Dynamic programming for pairwise overlap alignment • Problem definition • Input: sequences ࠵? (of length ࠵? ) and ࠵? (of length ࠵? ) • Output: the alignment between ࠵? and ࠵? involving subsequences of nucleotides that extend to either one end for both ࠵? and ࠵? in different directions, or both ends for either ࠵? or ࠵? whose alignment score is the highest • One end of a read aligns to the overhang of another (case #1 and #2 in the previous slide) • Part of one read aligns to the other entire read (case #3 and #4 in the previous slide)

Lastly, let’s look at overlap alignment between ࠵? and ࠵?. Similarities from global and local alignment: ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? !#$ , ࠵? "#$ + ࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵? • Initialization: the same as local alignment. Thus no penalties for indels in the beginning of the two sequences, i.e., 0’s for all initial entries ࠵? 0 , ࠵? = 1, … , ࠵?, ࠵? 1 , ࠵? = 1, … , ࠵? . • The back tracing starts at the end of either ࠵? or ࠵? , and ends at the start of either ࠵? or ࠵? . • For ࠵?(࠵? 0 , ࠵? 1 ) , it is the same as the global alignment. Thus Only this has never seen before.

Overlap alignment 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In overlap alignment, we use 0’s for all initializations. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? !#$ , ࠵? "#$ +࠵? ࠵? ! , ࠵? " ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Overlap alignment 0 0 0 0 0 0 0 0 0 0 0 ? ? ? 0 0 0 0 0 0 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In overlap alignment, we use 0’s for all initializations. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? !#$ , ࠵? "#$ +࠵? ࠵? ! , ࠵? " ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Overlap alignment 0 0 0 0 0 0 0 0 0 0 0 1 0 -1 0 0 0 0 0 0 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In overlap alignment, we use 0’s for all initializations. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? !#$ , ࠵? "#$ +࠵? ࠵? ! , ࠵? " ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵?

Now we try to fill in the whole matrix for overlap alignment.

Overlap alignment 0 0 0 0 0 0 0 0 0 0 0 1 0 -1 1 0 -1 -1 -1 -1 0 0 2 1 0 2 1 0 -1 -2 0 -1 1 1 0 1 3 2 1 0 0 -1 0 0 0 0 2 4 3 2 0 -1 -1 -1 -1 -1 1 3 5 4 0 -1 0 -1 -2 0 0 2 4 4 0 -1 -1 1 0 -1 -1 1 3 5 A C G A C T T T G A C T T T C G ࠵? ࠵? ࠵? ࠵? ’ , ࠵? " , ࠵? = 1 … ࠵? ࠵? ࠵? ! , ࠵? ’ , ࠵? = 1 … ࠵? Gap penalty = d In overlap alignment, we use 0’s for all initializations. We use ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? !#$ , ࠵? "#$ +࠵? ࠵? ! , ࠵? " ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵? Back tracing starts from the maximum entry in the rightmost column and bottom row in the matrix.

Overlap alignment Back tracing starts from the maximum entry in the rightmost column and bottom row in the matrix, and ends at the first row or first column in the matrix. In the corresponding alignment, we • don’t show indels that overhang the sequence; • do show indels within the alignment ACGACTTT-G ACTTTCG overhang

Four cases of back tracing in overlap alignment ࠵? ࠵? Start Rightmost column Last row Last row Rightmost column End First row First column First row First column 0 0 0 0 0 0 0 start end 0 0 0 0 0 0 0 start end 0 0 0 0 0 0 0 start end ࠵? ࠵? 0 0 0 0 0 0 0 start end Case #1 Case #2 Case #3 Case #4

Summary of overlap alignment • Initialization: • No penalty on the gaps in the beginning of the two sequences. • Thus ࠵? ࠵? ’ , ࠵? " = 0, ࠵? ࠵? ! , ࠵? ’ = 0 • Recursion or recursive function: • Back tracing: • Start from the entry/entries on the rightmost column or bottom row that has/have the highest score; • End at the first column or first row. ࠵? ࠵? ! , ࠵? " = max ࠵? ࠵? !#$ , ࠵? "#$ + ࠵?(࠵? ! , ࠵? " ) ࠵? ࠵? ! , ࠵? "#$ − ࠵? ࠵? ࠵? !#$ , ࠵? " − ࠵? like local alignment like global alignment unique for overlap alignment

Summary of pairwise alignment using DP. Global Local Overlap Initialization Linear indel penalties No penalty, all zeros No penalty, all zeros Sub-problems to problem Only diagonal, upper and left Additionally a zero in max Only diagonal, upper and left Start of back tracing The bottom right entry Any entry with the max. value Entry on the rightmost column or the bottom row that has the max. value End of back tracing The top left entry Whenever it hits a zero Whenever it hits the first row or the first column

Summary of dynamic programming algorithm • Definition: An algorithm to solve a big problem by breaking it into smaller sub-problems. • Application: only applicable to problems whose optimal solution of a big problem depends on the optimal solution of its sub-problems • For example, pairwise sequence alignment • Warning: Not every problem can be solved by dynamic programming, like maximum matching problem and finding the longest path problem in a graph. • What is essential? • Initialization (numbers on the first row/column of the matrix) • A process of constructing the optimal solution of the problem from its sub-problems • In our case, it is the equations. • Rule of back tracing

A little about computational complexity • It is important to know how many operations we have to compute for an algorithm. • It gives us some idea of how much time the algorithm will take to finish. • However, the number of operations depends not only the algorithm, but also on the inputs, i.e., in our cases, the sequences length, which are ࠵? and ࠵?. • To take into account the input size, we use big-O notations. For example, in dynamic programming algorithm for the three types of alignments we showed, each of them requires filling in a (࠵? +1)*(p+1) matrix, and each entry takes just one operation, which is a “max” calculation. So the big-O notations for these three alignments using dynamic programming algorithm is O(mp). • Here “+1” is eliminated as we usually consider the case when ࠵? or ࠵? is very large, and ”+1” is so small that they can be ignored compared with the majority of computations.

Which computational complexity is considered good? • Suppose ࠵? and ࠵? are both very large and that they are very close to each other, and here I use ࠵? to represent both of them. • The quadratic one, ࠵? ࠵? - , although can be tolerated, is not a very good algorithm. • In fact, anything bigger than ࠵? ࠵? - is horrible, as it will take months/years to run on a human genome, which is three billions base pairs long. • The linear one, ࠵? ࠵? is considered as very good. • The optimal one, ࠵? log(࠵?) is very hard to come up with, if not impossible. • The constant is almost nonexistent for sequence alignment.

Would you push the limit down and save computational time / lives?