Video 5-1: The read mapping problem in bioinformatics

0:02 So welcome back to uh algorithms of bionformmatics. Chapter 5 is about string matching. Biologically motivating is a slightly different term for roughly the same problem, the read mapping problem.

0:15 Uh and I'll uh I'll start commenting a bit on that and then we'll go through several rounds of formalizing the problem and solving it in different in different aspects. the same motivation will keep us uh busy for a bit longer than just this part. Uh because there's a few um um shortcomings of of the part that we'll we'll cover today that for biology aren't sufficient and also the performance wouldn't be good enough. So uh that's a it's a big chunk. Um also all developed in the last 25 years. So it's compared to some other algorithms that we we covered here and that you cover elsewhere is is fairly recent stuff.

1:04 What's the read mapping problem about? So you may remember from last uh unit that we were looking at genome assembly and uh the the exploding newspaper problem as a metaphor for what this puzzle is look is looking like. uh you have these many short fragments and you're trying to reassemble a genome from that. Uh the technical term for this approach by the way is shotgun sequencing and that's the official term that people use where the shotgun corresponds to the many small bits. You're somehow breaking your your DNA into these small pieces and you only read off the small small pieces. uh the puzzle the the burnt pieces of newspaper this the short segments are called reads in this context that's the pieces of DNA that you can read I think I've used that term in the previous uh unit as well uh but here it comes up in a slightly different context um but remember these are just these snippets of DNA in terms of the scale um if we stay with human genomes you'll know uh if you wake up at the middle of the night, the number of nucleotides in the human genome is roughly 3 billion.

2:21 It's a good number to keep in mind. Um to make any of these techniques from last time work at all, you probably want to cover every position 30 times. And people often do a bit more to be on the safe side because these I showed the videos last time about these sequencing machines, they often produce much more than that in just one go. So you might as well um but a reasonable lower estimate might be 30x coverage. So hovering every position should be covered by on average 30 little fragments. That means that's 300 million little reads each of a few hundred characters. So roughly if you don't compress it already 60 gigabyte of of raw data. That's just the reads. Now these you would want to feed into the algorithms from last week to assemble a genome. And indeed for genome sequencing also from uh from some species that have larger genomes than humans the main bottleneck is not the lab techniques.

3:28 The main bottleneck is the computational part. So luckily we don't always have to do that. We now have for humans at least a reference genome that we can somehow try to make use of. Whereas this was called denovo sequencing. So sequencing from scratch as if we knew nothing. Maybe we don't have to do that all the time. And there's another u motivation for why we may not want to do denovo sequencing all the time. That is how full genome sequences um sequencing is actually used for for med medical applications.

4:14 It turns out that having the full sequence of your genome is great but is also overwhelming for medical uh diagnosis. We wouldn't even know so well what to do with this full genome. What is much more helpful for very specific kind of diseases or genetic um mutations that are known to cause problems is very much more localized. So to make it more concrete um this specific example uh that's explained with a bit more context in the book is some uh form of of disability of mental retardation um called the odo syndrome and people trace down from from patients suffering from this and comparing their genomes that it's a single mutation that causes a specific protein to be short cut short. Essentially, we know now uh that this means there's a stop code on where there shouldn't be one. And these can be caused by changing a single letter in the genome. And this is indeed something that happens occasionally because a single mutation is not as unlikely as um larger changes. And biologists call it a single nucleotide polymorphism or a snip.

5:32 And these are in this case indeed causing or indicating this specific syndrome. Uh and there's other examples. So if there was a way for say a person to sequence their genome and just check for these specific single base replacements that often is is enough to detect uh certain genetic diseases. We don't really need the full uh sequencing of the genome. for many cases we wouldn't even make use of that.

6:10 So let's try not to do that and instead imagine the slightly simpler task. There's the long human reference genome and you're trying for your reads which this is not to scale. They're much shorter, right? You're just trying to find the best place where they match and maybe in this matching process you detect ah it matches but one character and then that can be by looking up where certain genes are located in the reference genome where we have all that information. We can try to identify what it check for specific uh genetic diseases. Most of the things that commercial providers do with your genome if you ask them to sequence it and then give you some some analysis of potential genetic problems that make may at some point how haunt you that's usually what they do and that's usually all they do. Uh that's my current state of knowledge anyways. So to formalize what this means for the computer science part for the uh uh computational question is the read mapping problem. We're given a text that's the reference genome one long sequence and we're suppose supposedly that's that's uh correct and done and then a whole collection of patterns short reads that's from the new patient that we're checking.

7:34 Now just to comment on the notation we talk about the same thing it's only come going to come up later again the text is of length n and the patterns are of length mr. So they could be different lengths. We typically think of these as being much much shorter than N. Remember N is three billion. M might be 200. But the number of reads, the little P, this might also be big. That was 200 million in the human ref in the human genome example.

8:07 uh I'll I'll write P as the vector as a as an array of all the patterns and each of the patterns itself is a is an vector or an array of num of U sorry of characters of nucleotides in our case. Uh so I have a to have a bit more convenient notation P R is one pattern and P is the collection of all patterns. Okay, so that's the input. Uh what are we trying to find? We're trying to find the best matches for these reads in the genome. Now what best match means is up for definition. This is a a template for different problems. The read mapping problem. It's not really one well-defined problem yet. Uh best match uh I guess the gold standard would be something based on our alignment uh algorithms in particularly these semi-local alignments. That just means you're allowed to slide it along uh the genome in any way. Um uh another another term was fitting alignment. That's what's uh used in the book. So you're trying to use all of the pattern, but you can zoom into a specific region of the genome.

9:14 Now that turns out to be very expensive and not really doable with the amount of data we have. we can try to be a little less um precise and say instead of an alignment that includes insertions, deletions, I'm only allowing mismatches. Then maybe that's still not quite doable. Um and instead we say if I'm if the best I can find for my pattern is a match with many mismatches still, then I'm also okay to not report that because it's it seems to be some some error. Maybe there was some uh contamination of the sequencing probe, but it's it's not actually part of the human genome that I'm interested. It's some some other dirt particle, whatever it was. Uh so I'm only interested in the matches that are close enough. And for the snips, that's exactly what we do. Uh single nucle nucleotide polymorphism would be a match of the read with a single mismatch. So this this seems very restricted but is is actually uh very useful for this application we just discussed.

10:22 The last one is least biologically correct but is computationally cheapest. That's only exact matches. And we'll focus on this last one for now for two reasons. First of all, it's easier and we have to start somewhere. But also these other variants can partially be built on this. In particular, the case with a single mismatch, you can maybe intuitively without going into much detail, if there's a long string that matches with a single mismatch, well, there's also going to be a large part of that string that matches without any mismatch. The mismatch is only one. It can't be everywhere.

11:05 But also uh that it simplifies the problem is of interest to us here. And so for the first part, forget about the read mapping problem in these other variants. It's always going to be about exact matches for now.

11:27 That problem is is generally known as string matching or pattern matching and uh is a wellstudied wellst studied problem. uh to be able to talk about it without making any ambitious statement uh ambiguous statements excuse me uh let's introduce some notation for that um I'll keep things general except for some cases where we'll restrict to the parameters that are interesting for our bionformmatics context all the strings come from some alphabet that's just the allowable letters and the size of the alphabet sometimes pops up that's this little sigma uh I think all the all the algorithms essentially assume that these are just numbers. So the the first n inte sigma integers of course for uh DNA that would be just the four uh nucleotides. It could also be the amino acids if you're talking about protein.

12:27 The methods are more generic and everything in a computer is a string. So uh uh often times um you can apply the same techniques. Uh but for for us you can think of sigma being really small and fixed. Uh then a string is denoted as the uh sigma to the n as a a tle of length n. Uh sigma star is this unspecified length. So union over all lengths. I don't know if we even need this. Uh that occurs epsilon is the empty string.

13:02 That's important. And then there's a bit of notation that essentially just says strings are arrays. And I use the same notation as for arrays. So uh this you've seen many times before. SN um starting at zero. Um sometimes we'll write the concatenation explicitly. I've you've seen that before. I think uh substrings is important. So talk about some uh starting index and an exclusive end index. And two specific terms that will pop up prefix and suffix. Prefix is a substring that starts at the beginning. Suffix is a substring that ends at the end. Maybe familiar to many. Um in case of doubt you can look at this slide.

13:53 Uh with that we can further narrow down the problem we're looking at. Um it differs from the read mapping problem specifically in that we're only having one pattern here. So it's the same same kind of text n characters and then we have one pattern um of length m in the read mapping context. We'll have many of these patterns. So for now we'll just assume you do one pattern after the other.

14:21 Okay. Then uh what you're looking for this this you can also get into debates about uh the easiest might be to think about finding the first match of the pattern in the text. Uh for read mapping probably you want to do all occurrences of the pattern in the text.

14:42 If the pattern does not occur at all then we return some dummy value. Right? Maybe also keep that in mind. Um the way the notation is defined, you can quickly say an a length m substring of the text at starting at i is is this. So if that's the same as the pattern, that's a match. Okay.

15:08 Now in in undergrad modules, you would spend a bit more time on explaining this one bullet, how to do it trivially. I'm not going to go there, but let me briefly say that you can you can solve the string matching problem by trying every starting position I that was would be uh allowed and then checking the next m characters do they match at this position with the pattern or not. Okay, so that that always works and it can find the first match or it can find any all matches uh as you like. The total number of character comparisons, so individual characters between string and pattern that we compare is roughly n times m. And that's going to be a bit too much for uh read mapping. Um for a single pattern, maybe still doable because m is not humongous. Um but we'll we'll not want to do this for the combined length of all the patterns.

16:07 That said, uh if you look at what some some standard libraries do, for example, the Java library, they essentially use this uh this brute force string matching method. Um Python's method is much more sophisticated and would maybe even be good enough for this bionformatics context. uh the Java one uh the one I've I've seen when I last checked anyways um would for the kind of data we have not work well because the typical thing that that uh is tricky for these string matching methods is if the alphabet is very small because then you have a lot of partial matches just for out of random luck and also long repetitions both within the pattern and in the in the genome. uh these tend to be potentially dangerous for the string matching methods. So even though um this this is the worst case running time for these simple methods, they often are better on average. But in the case of bionformmatics data, the worst case might actually happen or close to the worst case might actually happen.