Video 5-7: SNPs, alignment seeds

0:03 The last thing I want to look at again is uh this promise reduction from exact matching to having small number of mismatches. And in particular, we'll revisit these these snips, these single nucleotide polymorphisms.

0:21 Um but the the construction is indeed a bit more general. So in a way you could say this this kangaroo algorithm is fine if it wasn't for this overhead of pre-processing and indeed it looks simple from the way I presented it but you have to kind of take the the blackbox construction into account if you judge simplicity for very small K we can do another trick that doesn't have strong guarantees in the worst case but seems for uh bionformatics data to work reasonably well or very well. That's the following observation. Uh if you have a pattern and it's supposed to match some position in the text with at most k mismatches then uh let's draw this so it's immediately clear. That's my text. Uh the text continues but I'm I'm uh doing mismatches. So all that counts is some specific uh starting position in the text. Now let's say k is uh three. So there's some mismatches here.

1:30 I don't know where they are. But the statement says if that's all then there must be some reasonably large portion of the pattern and text that match exactly. And that's for example here all of these are equal. Maybe there's some more outside. What's the worst case for this? Well, if you balance out the mismatch parts precisely at uh regularly spaced intervals, then you get this bound. Or more precisely, let's draw it for a concrete example.

2:15 let's do the very extreme case first. That's how you can check the corner cases. Say if K is one and M is three. The worst that can happen is that this is the M mismatched character. But then I'm supposed to still have at least one character that I match. And if you insert that, that's what t comes out to 3 over2 rounded down. Uh in a larger example, say uh we pick characters like this. I could probably for the worst case append more things, but uh maybe we don't. Then let's count. So that's 369 and K is three. Now, the bound says I'm supposed to have at least 9 / 4 rounded down is is just two. And indeed, I have two characters that are matching. I could have more, right? If these are somehow more clumped together, the mismatched parts, the the red crosses, then I might have a longer stretch where the pattern and the text match exactly, but I'm guaranteed at least T characters, consecutive characters match. Clear? It's easy kind of pigeon hole argument. Uh these matches have to go somewhere and they get you get the shortest fragments if the mismatches are as evenly spread as possible. And that that gives you this one.

3:48 So we can use this observation by breaking our pattern into k + one substrings each of length t. And if there's something left over, you make the last one a bit longer. that doesn't matter. Uh all that matters is that if I'm finding exact occurrences for these uh fragments, there's K plus one of them. So if there's an occurrence where this pattern uh occurs with at most k mismatches, these k mismatches can only kill k out of these patterns or these sub patterns. There must be one left. and whatever whichever you pick it has length at least t. Okay, simple approach. Um and so you can use whatever algorithm from our section about exact matches to find these these uh smaller patterns and the occurrences where they occur. Let's call those the seeds of a potential match. And then you can try to extend this. Um, for example, if this is the portion of the text and you found that this fragment matches exactly, then you know that the pattern continues to the left and right a little bit.

5:05 So, you just compare the other characters and count how many mismatches you have. That's easy to do. Now the hope is that you don't have too many seeds that don't also correspond to some uh actual match you want to output. Uh and unfortunately we have no guarantee about that. Um you can easily come up with a bad instance at least a very degenerate one. You could say the text is just all A's and the pattern is almost that, but it starts with a bunch of non A's. Um, now say if we fix K to be, I don't know, two. So, we're allowing two mismatches. I'm going to be mean to not have any. I start with three B's and then I fill it up with A's. uh like this.

6:08 So, I'm going to break my pattern into these triples here. I'm going to find a ton of matches for these triple A's, but none of them can possibly lead to an actual match. So, the algorithm finds uh n * k matches or seeds, but they're all spurious. So in the worst case um so and and this is the number of matches you have to check each one so that the running time is n * m is back to the brute force version.

6:40 So this idea that I sold you as the the practical solution to this problem has a really really bad worst case uh but um it happens to not occur so much. uh I couldn't come up with a with an equally bad example that's not equally artificial say if you think about biological data there's all sorts of weird things happening but in terms of long repeats and so on uh but there's so much self similarity is rare and this seems to be what kills the algorithm here so the worst case is bad just as bad as the brute force method uh if you believe that your data is is more gentle to you, then these are rare and you can get by with much uh less running time. At least what we can say about this is that it's not missing any matches and it returns no wrong matches. It's still an exact or correct algorithm. It's just that the running time we don't know how to bound better than this uninteresting brute force bound.

7:49 But in in that sense um the uristic part is just about the speed up. Now we can also take the same idea finding seeds of very similar parts in this case exact matches and try to extend them in a different way. Uh this picture I had on the previous slide here that assumed I'm just checking I'm just counting how many mismatches I have. I could also start now building the dynamic programming matrix of a local alignment left and right. I'm still going to restrict the segment of the text I'm going to look at. Now say if you have um okay let let's let let me draw this properly uh so that it's not just this site remark here.

8:44 This is the entire text. And here is a snippet of my pattern that I found exactly matches the text at this part. I can try to extend this to say an edit distance at most k match.

9:12 Well, the K means I can have at most K gap symbols ever in my alignment between the pattern and the text. So the the pattern will continue in some form and the text continues a bit. In the worst case, I might have uh k gap symbols somewhere here. So I'm using from the text maybe Uh it's a bit hard to write because I didn't give this a name. Say this position is J prime. That's the starting position of this of this pattern snippet. J prime up to have a name for the length. It was T before like this.

10:02 Then this is covering J prime characters in P and we don't have to look more at more than J prime + K in the text. If I go further left, there's no way I can get rid of all those text characters within edit distance uh K. And same for the right. Um it's M minus J prime minus T + one minus plus k some something of that sort the number of characters to the right of that end point and plus k again. So I can take a a substring of the text that is essentially m plus 2k characters and then I can find the best alignment between t prime and p.

11:25 and the cost for that is essentially bounded by m and k. It doesn't involve n. So, uh the seed helped me to make this cost much much more local. And it it doesn't even matter all that much whether I'm extending it just with k mismatches or an alignment. It's a factor m, right? If I do it like this, the running time would be m ^2 + k.

12:03 So that's that's the good thing about this idea of having seeds for matches that you can extend them now also in a more general sense. Now unfortunately this uh this will find the best possible alignment locally but there might be alignments that are good overall. So good alignments between P and a substring of T that don't have any good enough seed just because they're they're uh globally good but they're enough spread out. Uh this is a bit harder to avoid with the general alignment scores anyways. So in uh in full generality this idea is uh then also uristic in the sense that it might miss some occurrences or uh if you if you tune it so that it finds all then uh it might spend more time. That said this is what many of the practical search tools are built on. um they often don't just do snips so uh mismatches but rather try to do local alignments. Um in particular there's two tools which are not targeted for read mapping but rather you found in some sample taken from some environment. you found some snippets of DNA or or a protein and you're trying to find what kind of uh known sequenced genomes of certain certain animals or or other living beings match this snippet that you found and for that a local alignment not a semi semi-local alignment is often more helpful but uh by and large they're trying to solve a similar problem you have a pattern you try to find that matching in a much much larger string or database of strings.

13:52 Um, and they use similar approaches um that have a a seed. Um, I'm I'm commenting on on the second part in in a bit. They have a seed and they extend this to an a local alignment between the two to the matching port of the text and the pattern.

14:12 You can even push the idea slightly further. Um, the seat doesn't have to be an exact match. But you can still use exact matching algorithms in the following way. If you restrict yourself to fairly short seeds and uh fairly tight thresholds for how much the seed is allowed to differ from the matching part of the text, then you can just premputee all possible strings, all ways to modify the pattern. Yep. so that uh this modified pattern and p are within distance sigma according to your alignment score whatever uh score you use that's in general an exponential set uh the longer p gets um I think k here was not a a good choice uh because this is let's say t k was the number of mismatches all the time um this this set of strings p prime will grow exponentially with t at least. So it's it's uh not feasible for large t but for small t you can do it and then you can search all those meddled width patterns exactly in the text. That's what this blast heristic is using and that's still the gold standard for people to search in in these uh databases of known sequences.

15:35 So they seem to have accepted the fact that this may may be a uristic in that it doesn't necessarily find all matches. Uh but it's the thing that is fast enough to scale to the size of data sets people have.

15:56 Okay. All that's left for this part for this uh unit is a quick summary what we've seen. Uh we spent some time on exact matching. Um and if you want to summarize all the different algorithms, you could say the running time is linear in n and m which could either be the length of one pattern or the total length of several patterns uh plus the size of the output which is only relevant for a whole chorus where the number of matches could be uh bigger than the text. You can have uh several patterns match at the same position. But this is essentially optimal that there's uh that's the best you can do in this setup. Uh whereas for all of the inexact matches, we only have partial results. For k mismatches, we can do the kangaroo algorithm in n times k time after this uh clever pre-processing step or we do these more uristic algorithms based on small seeds that match exactly. You want to go back to the original gold standard problem. You can either pay the price for the exact solution, usually not feasible, uh or you can from what we've seen so far use seeds and try to locally extend them into an alignment. But there we have very little formal guarantees.

17:12 We're not even sure we're finding the the the correct set of matches. And whether it's fast or not, you're you're relying on your assumption about the input. And I'm not aware of some some formal arguments in either direction what kind of inputs uh what what properties the input needs to have.