Video 7-8: String matching in BWT-compressed texts, backward search, FM Index

0:03 That's lovely. We did want to have uh a data structure for read mapping. If you remember from all all that time back, we started with the bionformmatics motivation of aligning small reads of DNA to the reference genome. We can't do that with that BWT. So that's not quite there yet. Turns out you can with a slightly different trick. And this is this is the most beautiful idea in in this part. So I hope we uh we can get it across in five minutes. It's actually an easy trick. You just search the pattern backwards and then you can use what the suffix array and BWT structure already provides you to uh narrow down the search in a very natural way. We could already do a binary search and that's all right. Um, but that only works uh if we have the suffix array. Now, with BWT alone, that doesn't seem to work. But we can do you could say a backwards radic search. Uh, fancy term for fancy name for something I haven't shown you. So, let's just see an example instead. Say we're we're looking for this small pattern Anna, which is a bad choice because it's a palendrome, but I wanted a small example.

1:23 uh initially we know nothing. So the range where we could be looking is all the suffixes. So all the positions in the BWT. What we really have is just the BWT for now. And uh maybe um that's that's all for the first step. Now suppose we start with the last A. We try to find occurrences of this pattern that at least have this last A. it has to be suffixes that have an a here and then we go further to the left. So the the end positions here we know are in that range in the sorted order because these are all the suffixes that start with an a. That's all occurrences of a full stop. That's not what we need in total, but for the last character that's fine. Now we include the second to last character and that means if this is the occurrence of a the previous character from there is wrap around the position in the BWT where there's an N. In this case there's only these two Ns. There could be several more, but I can take the first and last and shrink the interval that I'm searching in to everything that starts with a n because again in in the sorted order, these have to be adjacent.

2:48 Now note that there's three things that start with a n here. And if I take this n and look at uh where it where it points. So I need to find the second n and the third n in the occurrence. Sorry, sorry. Uh a n, not we need n a not a n. Uh there's this three ns down here, but only two that start with na. What I do is this is the second n in the bwt. it has to map to the second n in the sorted order. And same for this guy. This maps to the last n. So I know I can move my range if I look at n a as the pattern. I can move that from these four things to these two positions. And that's a range that contains everything that starts with na.

3:48 And I could do this from just again using the same rank and select queries on the BWT. In our case, the pattern goes on. So we have an A in front. So we check in the Bur Wheeler transform the range that we currently have highlighted. These are all A's. So I pick the first and the last A. These two in uh only. And that is the um first, second, the third and fourth A in the sorted order which corresponds to these two guys.

4:21 Okay. So here this maps to these two positions and that's that's my interval of occurrences of the pattern because I've uh reached the end of the pattern the beginning of the pattern now. So I know in sorted order I have occurrences of this pattern at ranks three and four. And I haven't done any binary search. I've just done steps through the pattern. And each step was fairly fast. Just this rangers queries log sigma. No no login. No n involved and anywhere. There's only one catch. I don't know where these occurrences in the text are. For that I need these two entries in the leaf array in the suffix array.

5:11 ignoring that part. Um, here's more carefully written code that does this backwards step. Um, it's doing the same thing. You initially start with the interval of everything that starts with the last character of the pattern. And then in each step, you look at the character of the pattern one position to the left and you find in that current range the first and last occurrence of that character in the BWT and map that back to the sorted column. Okay. And once you're done through all the pattern, you can return this interval of ranks. These are all the ranks where if your rank is in there, the entry in the suffix array would give you the starting position. So the only thing that's um missing is to do the same sampling trick for the suffix array. Previously we said we use the rank array and store every teeth. Um here we're doing essentially the same trick for the suffix array. And then by for each of these ranks in the interval decoding a bit more until we reach um an entry where we have the starting position stored. We can report where all these occurrences in the text are.

6:33 And the total time for this is m * log sigma. Again sigma tiny in our case. So it's essentially like uh in the suffix tree just uh m to find for example how many occurrences there are in the text and then here this this o is the number of number of occurrences. Uh if you want to find out where these are in the in the text the starting indices then for each of them you have to continue maybe for t more steps and read the corresponding value of the suffix array.

7:08 O okay that's it that thing is called the FM index and is at the root of widely used tools for doing read mapping. Um, if you have single genomes, it's still essentially state-of-the-art. Um, there's a bit of a gap between what I've shown you and the FM index and the read mapping because we still need to do inexact uh, searches. This is all for precise string matching, but we've talked a little bit about how to get from there to the inexact matches. Um the FM index was one of the first of these compressed self indexes and still widely used. Uh you can prove with the right wavelet tree representation that it compresses the text in something better than a Huffman coding trick. Um kith order empirical entropy for those who know. If not, forget it. Uh, and it's just the reference for the the read mapping tool that's based on it. Um, and this is an area that's still somewhat active. Uh, because the data keeps keeps growing and evolving. Even that trick is not good enough anymore for these 100,000 genome collections. You need to do um even even more aggressive compression.

8:28 And there's follow-up things on that, but uh we'll leave it with that. uh what's also not kind of a theory question might be can you do all the operations that you have in a suffix tree using just the BWT uh in particular we're missing the LCP array so far and for just the string matching you didn't need it for all the other fancy things that a suffix or enhanced suffix array can do uh we don't know yet if that's doable without extra space for these uh for these Self index indexes