Video 3-9: Multiple sequence alignment

0:03 after we've seen lots of different variants of pair wise alignments and ways to solve them um small polyilock speedups space efficient speed ups etc. At least the conclusion in all these cases was we have reasonably efficient algorithms for solving uh all variants of these problems. Now the um biologically more relevant version of the problem is not taking two sequences but taking many sequences that's called multiple sequence alignment MSA and just to briefly motivate again where this uh comes from in the biological context. uh one one one benefit or one one strength of pair wise alignments is actually a a problem in the applications that we use it for in that sometimes if you just take two sequences they will have some similar regions just because of random chance and then we would identify these as oh these are highly conserved these must be important etc when it is just a statistical fluke and this becomes much less likely if you don't use just two sequences as a kind of obvious obvious observation.

1:21 Um if something's just shared between two sequences like one sequence you fix the other is random. So you just have to match one random thing. If this happens simultaneously for many sequences uh is much less likely to occur in at least our completely random null model. So it has more significance and some bionformat bioinformatician put this in um much more metaphorical language.

1:48 Um here's an example from hemoglobin again that we've seen at the very beginning of this chapter. Remember there's two types of um of globin protein molecules that make up. So two copies of the alpha, two copies of the beta make one hemoglobin molecule in humans in healthy humans I guess. Uh if you look at this ver one version of these um in different species so there's us but there's also uh mice and chicken and zebra fish and this uh African clawed frog uh as some example of species. I mean, these are all animals, but they're not super close, right? They're not all mammals, etc. Still, you can build a multiple sequence alignment of the beta globin in all these species. They all need to transport oxygen. Even insects have hemoglobin, it seems. Uh, so this is a very widely shared function. And indeed all these um these proteins can be aligned with reasonable um agreement. So the star means exactly equal. No, it's not true. Is that it's a bit hard to see because I had to scale it down. I think the star should mean equal on all columns. It's just here.

3:13 This is the L, not the A. Uh so star means exactly the same in all five rows. um and colon and dot are varying degrees of uh still chemically similar.

3:30 So the colon still means chemically similar. I think the dot is essentially a mismatch and then there's also spaces for even worse mismatches or gaps.

3:41 The point is uh here you can see some some regions right here or uh I think it's surprisingly not not actually the most the longest stretch of uh conserved things you you can see that there are conserved regions And the hypothesis at least should be that these are these are important for the function. Whereas some other parts maybe here it seems less important to have exactly the right amino acids. And if you take any pair of the sequences you still get a kind of the same picture but it's not as clear. The signal for having this very strongly preserved region here is much stronger if you have more sequences. That's all um I want to say here as motivation. Um the biological uses of this are manifold. One is to take sequences where you already know they have a similar function and then you can zoom in where exactly is the interesting thing in the sequence happening. You can also take uh sequences where you where you know they should be similar but maybe uh there's variation. Um we didn't talk about ukareotic genes. So um let me skip that detour. Uh you can use it for for various things. You can also try to take for proteins where you know that they are highly conserved. You can take the similarities that you find to build um kind of a relation um a weighted relationship. Here you see that the sequences are actually grouped so that uh closer the more more similar sequences are closer together. And so it found out that uh human and and mice are closer whereas chicken is already a bit further apart and then there's a jump to zebra fish and these frogs for us. Uh what changes from the previous part in this chapter is we have k sequences instead of just two again over some alphabet. They may be of different length. We'll for the analysis often assume they're the same length just to keep notation simple. But you can definitely deal with whatever lengths an alignment. Now it used to be columns two uh two letters in in a column where you could have a gap symbol. Now there's k of these. So I'm taking um the alphabet add this gap character that's my sigma minus or sigma gap and uh I'm taking columns of k of these symbols and then a sequence of these columns. That's an alignment. Right? So in this picture any of these uh five topples here that would be yeah not very nicely drawn that would be one of these entries and then the alignment is just a string of these columns.

6:46 The other thing that changes is we need a score. We need a way to decide um how good a certain column is. And there's actually many options. Um we've already seen this in a restricted form in the last chapter where we talked about motives. There's this picture from they also align sequences on top of each other and we uh we said we have to somehow assign a score for how well this alignment is. uh the difference with motives there it was um this promoter region in front of a gene where remember maybe the uh the video over these uh angry looking blocks that you can't pass um no RNA polymerase for me um here it's really known that they are preserved in exactly the right length there might be single replacements of of nucleic acid bases uh but there's no insertions or deletions possible really in these motives for the general multiple sequence alignments bull as as before you can stretch things by adding these gap symbols.

7:48 So you could take um the score that we've used in in unit two the total hamming distance which as a recap you kind of find a consensus string by doing a majority vote in each column and then you count how many rows are different from the consensus. Uh here we'll do something different which maybe a bit better reflects um the scores if you have gaps in there and uh if you expect that maybe the sequences are not really very equal here for the motives the idea was there should be extremely conserved parts. Um here we're a bit more relaxed. So what we what we'll use is called the S spre the sum of pairs score. That's where the SP comes from. Uh so it's again you sum up the scores of all the columns and the score of one column with these K letters is just summing up all pairs in that column of alo all all possible row row pairs and then using any pair-wise alignment score. So the ba main benefit of this SP score is you can use whatever fancy pair wise alignment score you've worked out to be useful and you could lift it to a multiple sequence alignment.

9:08 Okay. So again the total score for one alignment you would have to sum this over all the columns. That's the same as before. It was always additive.

9:27 Let's first look at how we can compute those um in the same way generalizing what we've done so far. Uh this is one one picture from the book um where they illustrate again the 2D version pair wise alignment and I just want to stress uh you can think of it as finding a shortest path from topmost to bottommost but the edges have labels have weights in particular for edit distance the diagonal ones where there's a match like this red one here has a weight of zero they're for free whereas as um a diagonal edge like uh this one where it's a and c as a cost one and same for the vertical and horizontal edges. But uh with these edge costs, you can think of this as being finding a shortest path in this directed asyclic graph.

10:21 If you look at more than two sequences, you can do the very same thing. It just now becomes a kdimensional uh graph. So instead of uh this this rectangle it becomes cuboid or whatever k dimensional um thing you want to draw there for three I can still draw it I guess maybe it looks like this and then you would I don't know start up here end up here and you have to somehow make your way through there that's fine so the same kind of dynamic programming approach works Now if you go higher and higher in the dimension at each point you can go uh in the same stay in that dimension in the same or go one step back but you can do this in all combinations. So you have 2 to the k different combinations um minus the one where you don't do anything. So here we had three options which is 2 to the 2 minus one.

11:26 uh with three sequences we already have seven options 2 to the 3 minus one and this blows up exponentially clearly. Another way to say the same thing if you count the number of vertices in this graph that we're trying to find shortest paths in that number is is roughly n to the k if the strings have length n and we have k strings. So you can do this for a few uh or a few sequences. Maybe a a handful is fine, but quickly explodes.

11:58 So that seems like we can't do this. And maybe not so surprising if you think back to what happened with the motive finding uh problem. This is an NP hard problem in in the general case. Uh so specifically this um okay the the the the general version proven in this in this paper is whatever um whatever edit distance score you take so whatever pair-wise sequence score s that you take as long as it's a metric and uh whatever alphabet size you have as long as it's two letters then the problem to find the best possible alignment um where the number of strings is part of the input and the thing you optimize is the SP score. That is an NP hard problem. Um maybe I can morally sketch the reduction briefly. Uh but I will I'll refer you to the paper if you want to see the calculations. I think the reduction itself is is cute enough. Uh but the calculations get messy. I think uh the reduction is from vertex cover on cubic graphs. So cubic just means three regular. So every vertex has exactly three edges.

13:32 Uh so we're taking input for this problem a graph of this sort and a number and the question is can we find a vertex cover right and this is supposed to be a cubic graph so regular okay Now we're constructing some strings and some thresholds and uh we'll actually use edit distance.

14:14 So um the distance we we choose is edit distance. Then the SP score is the sums of edit distances. Um we can use sigma is two. I'll only show this for the special case. Okay. So um maybe this doesn't belong into the proof here. I'm I'm pretending that I'm only looking at the special case and uh you can show the same construction also works in more general but the calculation becomes even nastier.

14:45 We'll construct three kinds of strings. Um, one are strings that are a one for every vertex and then there's lots of zeros in between. And the exact number you have to work out, but it's the same number between each each uh pair of ones. I'll just copy this.

15:17 It roughly matches. Maybe it doesn't quite work. So these are template strings and we will not have a single one of these but many many copies. So B copies of these Then we'll also do uh another thing. So this is one one for each vertex v_ub_1 v_ub_2 and so on.

15:58 Uh then we also have b copies of a string that forces us to select a subset of vertices of size c. And this is a string that has just c1s. Okay. But in any good alignment, you'll have to put these in alignment with uh some one columns because again you have many copies of these. The copies of these first string, they will all be aligned with each other. Anything else just drives up the cost by so much. If you make big B big enough, it'll uh be above the threshold. Uh same for these. you all all of these B copies have to be aligned to the same vertices. Anything else the cost already blows up too much. So these select columns of vertices. And now you have to uh add um for each edge in uh in the graph or undirected, right? Let's just write it like this. For each undirected edge in the graph, you now have a string that again looks a bit more like this first one.

17:11 Um, but it only has ones for the two vertices that this edge uh connects. So you may have one that has zero here. I'll I'll copy this whole thing. Maybe it's cleaner.

17:36 So these are all zeros except for say these two which are ones. But now I'm I'm going to shift this thing around a bit. Um there's going to be extra zeros at the beginning and at the end.

18:03 Uh, and between them I'm going to delete a bunch of zeros roughly the same. So this means um Unless I want to pay at both ends for all these zeros, I'll have to shift this one inward. Oops. In any good alignment, this will rather look like this. Because if I leave this sticking out at the end, it's going to be too expensive to match those zeros to replace this with gaps in all the other strings. So, I'll have to align it so that only as much sticks out left and right as needed, which means I can only put one of the ones in agreement with some of these picking vector ones. And you do this for every edge. And then you can maybe see that you can pick a threshold so that if there's a vertex cover, every edge can be aligned so that at least one of the two ones is in agreement with all these B copies.

19:37 Again, think of B as close to infinity. Um very very big uh polomially big. Um, if you can at least bring one of the ones in agreement with one of these picked vertices, then you're not paying as much. If you have any edge where you cannot p pick either of the two ones to align with one of the picked vertices, you'll pay this cost of uh all all the pair-wise um distances with with this column will pay extra. And that signals there's no vertex cover.

20:40 So let's call this is this is T strings. This is P strings. These are edge strings. And then the game becomes counting exactly what the SP score of this alignment is, showing that any other alignment can only be worse and showing that you can you can separate the yes and no cases and uh is is a bit messy. So I'll spare us the details.

21:16 Now we still want these multiple alignments. Uh if it's NP hard, what can we do? Someone should do a lecture on this. Um one option is you give up optimality and uh that's um something where the SP scores are fairly well suited for. Um first observation is if you have the uh spore of some k sequences a1 up to a k you can always upper bound it by producing some alignment right that's uh pick pick some alignment that you constructs in some some way I'll show you more in a second that's an upper bound for the the distance you get from the optimal alignment and you can also lower bound this sum of pairs alignment by taking the actual pair-wise distance ignoring the rest. So in in um I'll I'll have notation for this um on the next slide. Uh whenever you have a multiple alignment, you can just focus on two rows, forget about the rest. That is a valid alignment of these two strings and it will have some score uh evaluated with the the two the pair-wise score. It may not be the optimal alignment if you just focus on these two strings because you were constrained by matching all the other strings at the same time but it is some score. So whatever you do here um the sum of scores here sums up the same it's the same sum right. Um but here I'm taking uh a different score. It's the score of this this alignment where I only look at the i and j row. And this may not be the optimal.

23:07 Okay, the notation we'll use m for this. Sorry, this is a bit python like syntax. Uh if you think about this um taking all all columns but then picking out rows i's and i and j. What what it's supposed to mean is if you have a specific multiple alignment, in particular one that gives you the best SP score, you can take out the rows that correspond to strings I and J and take their score. That score is lower bounded by the optimal alignment of these two strings in isolation. If you forget the rest of the other strings, that's why for this sum of pair scores, you get an upper bound from any alignment and a lower bound from the best pair-wise alignment. And so there's some hope that you can approximate this and um use this in a in a branch and bound algorithm for finding shortest paths. Again, remember this kd dimensional uh uh edit distance graph.

24:07 If you somehow have a way to prune off entire branches of the graph, then you're in good shape. And that's a standard approach. If you have uh a lower bound for the distance for some node to the target and that is worse than what you know for some other vertex as the entire distance then you can skip even looking at that vertex and for that we need precisely some some lower and upper bounds for uh our true distance that we can quickly compute. So the last thing I want to um show you is that we can indeed get an efficient approximation for this SP score and that's useful in isolation but you can in particular also use it as building block for this branch and bound algorithm which is still going to be an exponential time in the worst case but oftentimes substantially faster than this n to the k.

25:06 Okay. Um, if you think about how you could potentially approximate multiple alignments, you somehow I don't know, maybe there's other ideas, but one one intuitive idea is that you somehow add one string at a time and hope for the best. Somehow uh treat this as an incremental problem. Add one more string each step.

25:34 That's at least something to worth trying. And uh sometimes it works. For example, here are three strings. Um they're not at all uh chosen specifically to illustrate a point. One has A's and then T's. Then this has T's and then G's. And this has A's and then G's. Along these edges, there's the optimal pair wise alignment where here you want to bring the A's into uh alignment. Here the T's and here the G's. And it turns out here you can combine these pair-wise alignments in a smooth way into this multiple sequence alignment and you keep all the matched pairs um in the multiple alignment.

26:18 Now this does not always work. If you take this last string and reverse it, right? Uh then this alignment here now wants to align the middle G's. This alignment wants to align the A's behind the G's and this wants to align the T's before the G's. So in a way we've created a cycle, right? Uh here we would want A then T, then G. Uh but here we want the T's, then the G's, then the A's, and here we want the other direction. So that that's not going to nicely go together. In a way, you can say these pair-wise alignments are incompatible. There's no way to just glue them together and uh still get a valid alignment. They um they split up strings in an incompatible way.

27:23 The thing that goes goes wrong here is precisely that we looked at things in a cycle and if you don't do that everything's fine. In particular, you can uh you can think of arranging your sequences in some tree. So for three you could have done just a path and in general you can do a path but uh uh it works for any tree. Then you can compare compute the best pair-wise alignment along all of the edges.

27:56 And then you can show whatever alignment you have between A1 and A2. Let's do some concrete example to make this bit more concrete. And maybe it's a good idea to just use this bad example here.

28:27 So we have to decide on an order. Instead of going all the three edges, let's delete arbitrarily one of them. Say this one. Then whatever order we pick, say we first take these two. So we start with this a a t and then four gaps.

28:58 So this was say a1, this was a2. And now uh we've we've cleared this edge. We've kind of contracted this edge, right? By building the alignment. And now we add in a3. And so we have an alignment with a2. I should probably write down uh let's not be confusing more than necessary. Let's say this is a1. This is a2. Now I'm trying to add a3. I have an alignment with a2 which is uh right here.

29:39 Now we'll need to make some room, right? The alignment I have with um A3 and A2 would say I have a a T. I have four gaps and then I'm putting my string A3 below here.

30:01 Now, uh, for making this a complete alignment, there's, um, things left here, right? Um, sorry, um, yeah, missing these four. There are some other columns that are gaps in A2 that didn't exist in our alignment for A2 and A3. That's these four.

30:30 For these, we just add another round of gaps in A3 and we still rem retain a valid alignment. And what's more, as long as we decide to treat the alignment of two gaps as zero, uh, this didn't occur in a pair-wise alignment. But in the multiple alignment, it can happen that there's the first string sticking out with more G's. But the pair wise alignment between A2 and A3, if I just zoom in on these two now, there's these gap gap columns and they'll they'll we we just treat them as a zero. That's my convention, but it also is a sensible convention. Then the statement is whenever I do this procedure with a tree then uh the score of the alignment that is induced by this multiple alignment for say A2 and A3 is the same as this optimal score that I started with because all I did to make it fit was adding alignment columns where both are zero. So I'm adding zero right now. Unfortunately, this does not mean that the resulting multiple alignment as a multiple alignment is optimal because there could have been some better way of uh grouping things and maybe there's no way of building a tree to give us the best alignment. All of these things are possible, but at least we have an algorithm to compute a multiple sequence alignment um that starts from optimal pair wise alignments and uh we can still choose what tree we want to pick.

32:23 Uh for the quick analysis that will come um a bit of notation to make things concrete. uh capital M will be the multiple alignment. So think of this as a as this matrix of letters or gap symbols.

32:36 Uh again formally it's it's a sequence of columns and each column consists of K letters from this blown up alphabet the actual alphabet and the gap symbol. DSP is the sum of pair scores and that's uh that's something to evaluate a multiple alignment and it depends on the pair wise alignment score that we'll use. Uh all of all of this you can just replace s with edit distance. Uh but the the ideas and algorithms are are more general. ds uh I think I've used it consistently like this before but that's the score of the optimal pair wise alignment between two sequences whereas s itself was a function that gives you the score of columns or the score of some alignment and you just give the alignment and it may or may not be the best alignment whereas ds d for distance should give you the score of the best alignment and then here's what I mentioned mentioned earlier if you have some multiple alignment you can just take two of the rows and ignore the rest that is a valid alignment for a i and aj if you ignore these double gapped columns and the score of these two rows will never be better than the optimal alignment but at least it's some alignment okay now to turn this into an algorithm we need to pick a And what's the maybe simplest possible tree? A star. Maybe you should have said a path, but uh uh for the analysis, it's a star. So there's some sequence that we put in the middle. A C for center and all the others are just around that one connected with some edge.

34:36 The algorithm starts by um computing all pair-wise distances because we use them anyways in the next step. So for every I and J we compute the optimal edit distance say then we find one of the sequences by essentially trying them all one after the other. Uh we compute what is the sum of the distances from that sequence to all the others. That's this sum. And we there will be one that is best. Let that be our center of the star. That's in a way the closest sequence to the others.

35:12 And then we can construct our multiple alignment following the previous algorithm by essentially merging in one uh one one leaf of the star at a time. I think the algorithm is fairly simple. All you're doing is compare pairwise alignments. This merging in really consists of just ignoring columns where there were additional gaps that didn't exist in your alignment. Right? This is all that happens. So it's very simple. But how good is this? Uh it turns out it's a two approximation or at least uh it's actually a slightly bit better when k is is small. You can show it's um 2 - 2 over k. Um let me sketch at least um how do you get this? And uh probably good to keep the notation somewhere visible.

36:23 There's one assumption here and that's also kind of a hint for the proof if you want. We need our pairwise uh distance that we base this on to be a metric. And I'm sure you're all uh eagerly working on the next assignment or finish where you'll show that unless you do stupid things uh global alignment distances uh tend to be um metrics. So this is not a strong re um requirement. All the sensible ones uh fulfill this. In particular, edit distance is a valid metric.

37:07 And then we'll show um the following. So we have some optimal alignment multiple sequence alignment and we have the thing that our algorithm computes the center star approximation.

37:34 And to make things uh to make the calculations nicer instead of sum of pairs I'm doing sum of ordered pairs. What do I mean by that? Um we said uh this SP distance corresponds to summing ordered pairs. And let's use V to say I'm summing Well, you can write it I comma J I guess all the distances. So here um what do I want? Uh it's the it's the score of the multiple alignment zoomed into I and J. Uh let's let's just say um in general these are for the two sums. zero. Um, yeah. No, no, no. Let's leave it.

38:50 uh v is a bit nicer to calculate but really if you're summing over all ordered pairs then you're just summing things twice and the entries where i and j are the same they'll contribute zero because they have distance zero okay so whatever you do b of m is going to be 2 * the sum of pairs distance of the same alignment in particular for an approximation guarantee where we take the ratio of MC divided by the cost for optimum that ratio is the same whether I multiply both things by two right so uh that doesn't make a difference okay uh first step is let's get a lower bound for the optimum so V of M star that's again it's the double sum i inj J the score of the alignment I get when I zoom in on rows I and J. Now this is some pair wise alignment of rows I of rows I and J. That's going to be never better than doing the pairwise alignment ignoring what the uh multiple alignment says.

40:16 That's going to be an a lower bound for opt. Then the second step will be to show an upper bound for our algorithm the cost for this MC and we'll basically show that it's let's call this um what can we call this let's call it little m I think that's still free uh we'll claim that this is essentially this approximation 2 - 2 over k * this m and then the result immediately follows from that by taking taking the ratio.

41:10 So why why would that be the case? what we get in the uh center star alignment. But first of all, same double sum, right?

41:30 That should be a star MC. I'm taking my column in and J. Now we'll have to go back to the um observation that I've mentioned here in the drawing. We didn't formally prove it. I mentioned it while pointing to this picture. Uh when you build an alignment from a tree like this, whenever you look at the pairs um that you had an edge for, their pair wise alignment is optimal. the pair wise alignment induced by the multiple sequence alignment only has some zero columns added otherwise its score is the same and because they're zero its score is the same as the optimal alignment.

42:24 So that's kind of fact from above or for all I the score if I take in my center star alignment a column C and this J.

42:49 The score of that thing is the actual best possible sequence alignment of the center sequence and AJ. Okay, that's that. Now the second thing is our distance is a metric and that also means it's a metric on these um pair wise alignments.

43:12 So let's draw the star again. If you look at these sans, there's two types of seamans. There's ones where one of the iO is is the C. So I'm I'm summing across one of the edges. And then there's the other ones where you're summing over two ed uh edges uh sorry uh two leaves say these two. Now we're we would be interested in the edit distance between those two in the direct way.

43:56 And this we don't know. But what we do know is that taking the detour via C. That's never going to make it shorter, right?

44:49 So for all the sumans I have two options. I'm either summing over an edge or I'm summing over one of these. Then I can upper bound it by summing over the two edges.

45:00 So how often do we sum over each edge? Uh we sum once for i equal the leaf, j equals the center. I guess another time for i equals the center and j equals the leaf. And then for every other leaf we sum it once.

45:39 Ah we sum it twice because it's both pairs. Okay. So that sum is is indeed every every edge occurs exactly 2 * k minus one times there's k sequences in total which means k minus two leaves and I'm summing it twice for every other vertex. Exactly. Uh it doesn't make a difference if it's the center of or another leaf. I made life more complicated. So we're summing over edge every edge uh 2k minus one times. And then what we sum are the edges from i to the edge here.

46:37 C I and that one is upper bounded or is actually equal um by the fact to the true distance we have to compare m which is summing all of these guys to this thing.

47:27 I think I should have expressed This also in terms of the sum of the distance to the centers. this is supposed to be K times. That's um And that's because um we chose the center to be uh the one where this sum is minimized. So uh I can uh I can pick any I can hold any J constant and sum over all the others.

49:24 I'm picking one of the J's as my as my potential center. Then what what remains of the sum is uh summing up the distances to all the other things. I'm summing up all the distances holding one fixed distances to all the other uh elements and uh the center is the one where this sum is minimal. So if you sum over all the other things that's in the that's never never smaller than k times the minimum you have k such sums. So you can upper bound this one as k * our small m then it matches with um almost the scaling was wrong. uh this one and then if you divide you take the ratio you get the approximation ratio. Okay.

50:24 Yeah. This this step um is important because that is the only place where we use that the center is actually the best center that we could have chosen. Other than that, all we needed was the triangle inequality and um bit of of rewriting.