Lecture 5: Library Complexity and Short Read Alignment (Mapping)

• Download English-US transcript (PDF)

The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high quality educational resources for free. To make a donation or view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.

PROFESSOR: OK, so welcome back to computational systems biology, I'm David Gifford. I'm delighted to be with you here here today. And today we're going to be talking about a topic that is central to modern high throughput biology, which is understanding how to do short read alignment, sometimes called read mapping. Now it's very important to me that you understand what I'm about to say today, and so I'm hopeful that you'll be uninhibited to raise your hand and ask questions about the fine points in today's lecture if you have any, because I'd be totally delighted to answer any questions and we have enough time today that we can spend time looking at one aspect of this problem and understand it thoroughly.

An associated topic is the question library complexity. How many people have heard of sequencing libraries before? Let's see a show of hands. OK, how many people have heard of read alignment before, read mapping? OK, great, fantastic. Let's start with what we're going to be talking about today. We're going to first begin talking about what a sequencing library is and what we mean by library complexity.

We'll then turn to what has been called a full text minute size index, sometimes called a burrows Wheeler transform index, a BWT index, an FM index, but this is at the center of most modern computational biology algorithms for processing high throughput sequencing data. And then we'll turn how to use that type of index for read alignment. So let's start now with what a sequencing library Is. Let's just say that you have a DNA sample, we'll be talking about various ways of producing said samples throughout the term.

But we're going to assume that we have a bunch of different DNA molecules. And I'll illustrate the different molecules here in different colors. And we have three different types of molecules here. Some molecules are duplicated, because as you know, typically, we're preparing DNA from an experiment where there are many cells and we can get copies of DNA from those cells or the DNA could be amplified using PCR or some other technique. So we have this collection of molecules, and to make a library, we're going to process it.

And one of the things that we'll do when we process the library is we'll put sequencing adapters on. These are short DNA sequences that we put on to the end of the molecules to enable them to have defined sequences at the ends which permits sequencing. Now, if somebody hands you a tube of DNA like this, there are a couple questions you could ask. You could check the DNA concentration to find out how DNA is there, you could run a gel to look at the size of the fragments that you're sequencing. We'll be returning to that later, but these are typically called the insert sizes of the library that you're sequencing, the total length of the DNA excluding the adapters.

But we could also ask questions about how complex this library is, because it's possible to run experiments where you produce libraries that are not very complex, where they don't have very many different types of molecules. Now that typically is a failure of the experiment. So an important part of quality control is characterizing the library complexity where we want to figure out here complexity is equal to 3. There are three different types of molecules. And we sample these molecules.

And when we sample them, we get a bunch of DNA sequence reads. And typically the number of reads that we get is larger than the complexity of the library. Here we have a total of 12 different reads. And when we sequence a library, we're sampling from it. And so the probability that we get any one particular molecule is going to be roughly speaking equal to 1 over c, which is the complexity. And thus, we could use the binomial distribution to figure out the likelihood that we had exactly four of these type one molecules.

However, as n the number of sequencing reads grows to be very large, typical numbers are a hundred million different reads, the binomial becomes cumbersome to work with. And so we typically are going to characterize this kind of selection process with a different kind of distribution. So one idea is to use a Poisson, where we say that the rate of sequencing is going to be n over c. And we can see that here shown on the slide above is the same process where we have the ligation of the adapters.

We have a library and we have reads coming from the library. We have a characterized library complexity here, there are four different types of molecules. And the modeling approach is that assuming that we have c different unique molecules, the probability that we'll get any one of them when we're doing the sequencing is 1 over c. And if we do end sequencing reads, we can find out the probability that we'll get a certain number of each type of molecule. Let's just stick with the first one to start, OK?

Now part of the challenge in analyzing sequencing data is that you don't see what you don't sequence. So things that actually occur 0 times in your sequencing data still may be present in the library. And what we would like to do is from the observed sequencing data, estimate the library complexity. So we have all of the sequencing data, we just don't know how many different molecules there are over here. So one way to do with this is to say that let us suppose that we make a histogram of the number of times we see distinct molecules and we're going to say that we can observe molecules that are sequenced or appear l times up through r times.

So we actually can create a version of the distribution that characterizes just a part of what we're seeing. So if we do this, we can build a Poisson model and we can estimate lambda from what we can observe. We don't get to observe things we don't see. So for sure, we know we can't observe the things that are sequenced 0 times. But for the things that are sequenced at least one time, we can build an estimate of lambda. And from that estimate of lambda, we can build an estimate of C.

So one way to look at this is that if we look at the total number of unique molecules that we sequence, which is equal to m, then the probability that we observe between l and r occurrences of a given individual sequence times c is going to be equal to the total number of unique molecules that we observe. Another way to look at this is the very bottom equation where we note that if we look at the total complexity of the library and we multiply it by 1 minus the probability that we don't observe certain molecules, that will give an estimate of the total number of unique molecules that we do see.

And thus we can manipulate that to come up with an estimate of the complexity. Are there any questions about the details of this so far? OK, so this is a very simple model for estimating the complexity of a library based upon looking at the distribution of reads that we actually observe for quality control purposes. And let us suppose that we apply this to thousands genomes data, which is public data on human. And suppose we want to test whether this model works or not, so what we're going to do is we're going to estimate the library complexity from 10 percent of the sequencing reads, so we'll pick 10 percent of the reads of an individual at random, we'll estimate the complexity of the library, and then we'll also take all of the region the individual and estimate the complexity.

And if our estimator is pretty good, we should get about the same number, from 10 percent of the reads and from all of the reads. Will people go along with that? Think that seems reasonable? OK, so we do that. And this is what we get. And it's hard to see the diagonal line here, but there's a big oops here. And the big oops is that if we estimate the library complexity from just 10 percent of the reads, it's grossly underestimating the number of unique molecules we actually have. In fact, it's off by typically a factor of two or more.

So for some reason, even though we're examining millions of reads in this subsample, we're not getting a good estimate of the complexity of the library. Does anybody have any idea what could be going wrong here? Why is it that this very simple model that is attempting to estimate how many different molecules we have here based upon what we observe is broken? Any ideas at all? And please say your name first.

AUDIENCE: I'm Chris.

PROFESSOR: Hi Chris.

AUDIENCE: Is it because repeated sequences, so there could be a short sequence at the end of one molecule that's the beginning of another one, middle of another one, so [INAUDIBLE].

PROFESSOR: Chris, you're on the right track, OK? Because what we have assumed at the outset was that all of these molecules occur with equal probability. Right? What would happen if in fact there are four copies of this purple one and only two copies of the other molecules? Then the probability of sampling this one is going to be twice as high as the probability of sampling one of these. If there's non uniformity in the original population, that's going to mess up our model big time.

And that could happen from repeated sequences or other kinds of duplicated things, or it could be that there's unequal amplification. It might be that PCR really loves a particular molecule, right, and amplifies that one a lot, and doesn't amplify another one that's difficult to amplify. So somewhere in our experimental protocol pipeline, it could be that there's non uniformity and thus we're getting a skewness to our distribution here in our library. So the other thing that's true is in a Poisson, lambda, which is equal to the mean, is also equal to the variance.

And so our Poisson's only one knob we could turn to fit the distribution. So coming back to this, we talked about the idea that the library complexity still may be four but then there may be different numbers of molecules of each type. And here's an idea for you, right? The idea is this. Imagine that the top distributions are the number of each type of molecule that are present. And it might be that our original assumption was that it was like the very top, that typically there are two copies of each molecule in the original sequencing library, and that's a fairly tight distribution.

But it could be, in fact, that the number of molecules of each type is very dispersed. And so if we look at each one of those plots at the top, the first four, those are going to be our guesses about the distribution of the number of copies of a molecule in the original library. And we don't know what that is, right? That's something we can't directly observe, but imagine that we took that distribution and used it for lambda in our Poisson distribution down below for sampling.

So we have one distribution over the number of each type of molecule we have and we have the Poisson for sampling from that, and we put those two together. And when we do that, we have the Poisson distribution at the top, the gamma distribution is what we'll use for representing the number of different species over here and their relative copy number. And when we actually put those together as shown, we wind up with what's called the negative binomial distribution, which is a more flexible distribution, it has two parameters.

And that negative binomial distribution can be used, once again, to estimate our library complexity. And when we do so, we have lamba be the same, but k is a new parameter. It measures sort of the variance or dispersion of this original sequencing library. And then when we fit this negative binomial distribution to that 1,000 genomes data, it's going to be hopefully better. Let's start with a smaller example. If we have a library that's artificial with a known million unique molecules and we subsample, it gives you 100,000 reads, you can see that with different dispersions here in the left, k with different values from 0.1 to 20, the Poisson begins to grossly underestimate the complexity of the library as the dispersion gets larger, whereas the negative binomial, otherwise known as the GP or gamma Poisson, does a much better job.

And furthermore, when we look at this, in the context of the thousand genomes data, you can see when we fit this how much better we are doing. Almost all those points are almost exactly on the line, which means you can take a small sampling run and figure out from that sampling run how complex your library is. And that allows us to tell something very important, which is what is the marginal value of extra sequencing. So for example, if somebody comes to you and they say, well, I ran my experiment and all I could afford was 50 million reads. Do you think I should sequence more? Is there more information in my experimental DNA preparation?

It's easy to tell now, right? Because you can actually analyze the distribution of the reads that they got and you can go back and you could estimate the marginal value of additional sequencing. And the way you do that is you go back to the distribution that you fit this negative binomial and ask if you have r more reads, how many more unique molecules are you going to get? And the answer is that you can see that if you imagine that this is artificial data, but if you imagine that you had a complexity of 10 to the 6 molecules, the number sequencing regions is on the x-axis, the number of observed distinct molecules is on the y-axis, and as you increase the sequencing depth, you get more and more back to the library.

However, the important thing to note is that the more skewed the library is, the less benefit you get, right? So if you look at the various values of k, as k gets larger, the sort of the skewness of the library increases, and you can see that you get fewer unique molecules as you increase the sequencing depth. Now I mention this to you because it's important to think in a principled way about analyzing sequencing data. If somebody drops 200 million reads on your desk and says, can you help me with these, it's good to start with some fundamental questions, like just how complex is the original library and you think that these data are really good or not, OK?

Furthermore, this is a introduction to the idea that certain kinds of very simplistic models, like Poisson models of sequencing data can be wrong because they're not adequately taking into account the over dispersion of the original sequencing count data. OK, so that's all there is about library complexity. Let's move on now to questions of how to deal with these reads once we have them. So the fundamental challenge is this.

I hand you a genome like human. 3 times 10 to the ninth bases. This will be in fast a format, let's say. I had you reads. And this will be-- we'll have, say, 200 base pairs times 2 times 10 to the eighth different reads. And this will be in fast q format. The q means that there are-- it's like fast a except that our quality score's associated with each particular base position. And the PHRED score which is typically used for these sorts of qualities, is minus p minus 10 times log base 10 of the probability of an error.

Right, so a PHRED score of 10 means that there's a 1 in 10 chance that the bases is an error, a PHRED score of 20 means it's one in a 100, and so forth. And then the goal is today if I give you these data on a hard drive, your job would be to produce a SAM file, a Sequence Alignment and Mapping file, which tells us where all these reads map in the genome. And more pictorially, the idea is that there are many different reasons why we want to do this mapping. So one might be to do genotyping. You and I differ in our genomes by about one base in a thousand. So if I sequence your genome and I map it back or align it to the human brain reference genome, I'm going to find differences between your genome and the human reference genome.

And you can see how this is done at the very top where we have the aligned reads and there's a G, let's say, in the sample DNA, and there's a C in the reference. But in order to figure out where the differences are, we have to take those short reads and align them to the genome. Another kind of experimental protocol uses DNA fragments that are representative of some kind of biological process. So here the DNA being produced are mapped back to the genome to look for areas of enrichment or what are sometimes called peaks.

And there we want to actually do exactly the same process, but the post processing once the alignment is complete is different. So both of these share the goal of taking hundreds of millions of short reads and aligning them to a very large genome. And you heard about Smith Waterman from Professor Berg, and as you can tell, that really isn't going to work, because its time complexity is not going to be admissible for hundreds of millions of reads.

So we need to come up with a different way of approaching this problem. So finding this alignment is really a performance bottleneck for many computational biology problems today. And we have to talk a little bit about what we mean by a good alignment, because we're going to assume, of course, fewer mismatches are better. And we're going to try and align to high quality bases as opposed to low quality bases and note that all we have in our input data are quality scores for the reads.

So we begin with an assumption that the genome is the truth and when we are aligning, we are going to be more permissive of mismatches in read locations that have higher likelihood of being wrong. So is everybody OK with the set up so far? You understand what the problem is? Yes, all the way in the back row, my back row consultants, you're good on that? See, the back row is always the people I call on for consulting advice, right? So yeah. You're all good back there? Good, I like that, good, that's good, I like that, OK.

All right. So now I'm going to talk to you about what are the most amazing transforms I have seen. It's called the Burrows Wheeler Transform. And it is a transform that we will do to the original genome that allows us to do this look up very, very quickly. And it's worth understanding. So here's the basic idea behind the Burrows Wheeler Transform. We take the original string that we want to use as our target that we're going to look things up in, OK, so this is going to be the dictionary looking things up in and it's going to be the genome sequence.

And you can see the sequence on the left hand side, ACA, ACG, and the dollar sign represents the end of string terminator. OK. Now here's what we're going to do. We take all possible rotations of this string, OK? And we're going to sort them. And the result of sorting all the rotations is shown in the next block of characters. And you can see that the end of string character has the shortest sorting order, followed by A, C, and G and that all the strings are ordered lexically by all of their letters.

So once again, we take the original input string, we do all of the possible rotations of it, and then we sort them and wind up with this Burrows Wheeler Matrix as it's called in this slide, OK? And we take the last column of that matrix and that is the Burrows Wheeler Transform. Now yo might say, what on earth is going on here? Why would you want to take a string or even an entire genome? We actually do this on entire genomes, OK? Consider all the rotations of it, sort them, and then take the last column of that matrix. What could that be doing, OK?

Here's a bit of intuition for you. The intuition is that that Burrows Wheeler Matrix is representing all of the suffixes of t. OK, so all the red things are suffixes of t in the matrix. And when we are going to be matching a read, we're going to be matched it from its end going towards the beginning of it, so we'll be matching suffixes of it. And I'm going to show you a very neat way of using this transform to do matching very efficiently. But before I do that, I want you to observe that it's not complicated.

All we do is we take all the possible rotations and we sort them and we come up with this transform. Yes.

AUDIENCE: What are you sorting them based on?

PROFESSOR: OK, what was your name again?

AUDIENCE: I'm Samona.

PROFESSOR: Samona. What are we sorting them based upon? We're just sorting them alphabetically.

AUDIENCE: OK.

PROFESSOR: So you can see that if dollar sign is the lowest alphabetical character, that if you consider each one a word, that they're sorted alphabetically, OK? So we have seven characters in each row and we sort them alphabetically. Or a lexically. Good question. Any other questions like that? This is a great time to ask questions, because what's going to happen is that in about the next three minutes if you lose your attention span of about 10 seconds, you're going to look up and you'll say, what just happened? Yes.

AUDIENCE: Could you explain the suffixes of t?

PROFESSOR: The suffixes of t? Sure. Let's talk about the suffixes of tr. They're all of the things that end t. So a suffix of t would be G, or CG, or ACG, or AACG, or CAACG, or the entire string t. Those are all of the endings of t. And if you look over on the right, you can see all of those suffixes in red. So one way to think about this is that it's sorting all of the suffixes of t in that matrix. Because the rotations are exposing the suffixes, right, is what's going on. Does that make sense to you? Now keep me honest here in a minute, OK, you'll help me out? Yes. Your name first?

AUDIENCE: [INAUDIBLE]. [INAUDIBLE].

PROFESSOR: [INAUDIBLE].

AUDIENCE: What is dollar sign?

PROFESSOR: Dollar sign is the end of string character which has the lowest lexical sorting order. So it's marking the end of t. That's how we know that we're at the end of t. Good question. Yes.

AUDIENCE: Can you sort them non-alphabetically, just different ways to sort them [INAUDIBLE] algorithm?

PROFESSOR: The question is, can you sort them non alphabetically. You can sort them any way as long as it's consistent, OK. But let's stick with alphabetical lexical order today. It's really simple and it's all you need. Yes. in red is the suffixes in the last colored group on the right?

PROFESSOR: No, no.

AUDIENCE: What's in red?

PROFESSOR: What's in red are all the suffixes of T on the very far left. OK?

AUDIENCE: On the right, last column group?

PROFESSOR: The right last column group. That last column in red, that is the Burrows-Wheeler Transform, read from top to bottom. OK? And I know you're looking at that and saying, how could that possibly be useful? We've taken our genome. We've shifted it all around. We've sorted it, we take this last thing. It looks like junk to me, right?

But you're going to find out that all of the information in the genome is contained in that last string in a very handy way. Hard to believe but true. Hard to believe but true. Yes. Prepare to be amazed, all right?

These are all great questions. Any other questions of this sort? . OK. So, I'm going to make a very important observation here that is going to be crucial for your understanding. So I have reproduced the matrix down on this blackboard. What? That's usually there under that board, you know that. You guys haven't checked this classroom before, have you? No. It's always there. It's such a handy transform.

So this is the same matrix as the matrix you see on the right. I'm going to make a very, very important assertion right now. OK? The very important assertion is that if you consider that this is the first a in the last column that is the same textual occurrence in the string as the first a in the first column. And this is the second a in the last column, that's the same as the second a in the first column. And you're going to say, what does he mean by that?

OK, do the following thought experiment. Look at the matrix. OK? And in your mind, shift it left and put all the characters on the right hand side. OK? When you do that, what will happen is that these things will be used to sort the occurrences a on the right hand side.

Once again, if you shift this whole thing left and these pop over to the right, then the occurrence of these a's will be sorted by these rows from here over. But these are alphabetical. And therefore they're going to certain alphabetical order. And therefore these a's will sort in the same order here as they are over there.

So that means that when we do this rotation, that this textual occurrence of a will have the same rank in the first column and in the last column. And you can see I've annotated the various bases here with their ranks. This is the first g, the first c, the first end of line, end of string character. First a, second a, third a, second c. And correspondingly I have the same annotations over here and thus the third a here is the same lexical occurrence as the third a on the left in the string, same text occurrence.

Now I'm going to let that sink in for a second, and then when somebody asks a question, I'm going to explain it again because it's a little bit counterintuitive. But the very important thing is if we think about textual recurrences of characters in that string t and we put them in this framework, that the rank allows us to identify identical textual recurrences of a character.

Would somebody like to ask a question? Yes. Say your name and the question, please.

AUDIENCE: Dan.

PROFESSOR: Dan.

AUDIENCE: So in your original string though those won't correspond to the same order in the transformed string. So like the a's in the original string in their order, they don't correspond numerically to the transformed string.

PROFESSOR: That's correct. Is that OK? The comment was that the order in BWT, the transform is not the same as the order in the original string. And all I'm saying is that in this particular matrix form, that the order on the last column is the same as the order in the first column for a particular character. And furthermore, that these are matching textual occurrences, right?

Now if I look at a2 here, we know that c comes after it, then a, then a, and c and g, right? Right. OK, so did that answer your question that they're not exactly the same?

AUDIENCE: Yes. I don't understand how they're useful yet.

PROFESSOR: You don't understand how it's useful yet. OK. Well, maybe we better get to the useful part and then you can-- OK. So let us suppose that we want to, from this, reconstruct the original string. Does anybody have any ideas about how to do that? OK.

Let me ask a different question. If we look at this g1, right? And then this is the same textual occurrence, right? And we know that this g1 comes right before the end of character, in end of string terminator, right? So if we look at the first row, we always know what the last character was in the original string. The last character is g1, right? Fair enough? OK

Where does g1 would occur over here? Right over here, right? What's the character before g1? c2. where is c2 over here? What's the character before c2? a3. What's the character before a3? a1. Uh, oh.

Let me just cheat a little bit here. a1 a3 c2 g1 $. So we're at a1, right? What's the character before a1? c1, right? What's the character before c1? a2. And what's the character before a2? That's the end of string. Is that the original string that we had? Magic. OK? Yes.

AUDIENCE: Wouldn't it be simpler to look at-- to just remove the dollar sign [INAUDIBLE] or do you mean reconstruct from only the transformed?

PROFESSOR: We're only using this. This is all we have. Because I actually didn't use any of these characters. I was only doing the matching so we would go to the right row. Right? I didn't use any of this. And so, but do people understand what's going on here? If anybody has any questions, now is a great time to raise your hand and say-- here we go. We have a customer. Say your name and the question, please.

AUDIENCE: My name is Eric.

PROFESSOR: Thanks, Eric.

AUDIENCE: Can you illustrate how you would do this without using any of the elements to the left of the box?

PROFESSOR: Absolutely, Eric. I'm so glad you asked that question. That's the next thing we're going to talk about. OK, but before I get to there, I want to make sure, are you comfortable doing it with all the stuff on the left hand side? You're happy about that? OK. if anyone was unhappy about that, now would be the time to say, I'm unhappy, help me. How about the details? Everybody's happy? Yes.

AUDIENCE: So, you have your original string in the first place, though, so why do you want to create another string of the same length? Like, how does this help you match your read?

PROFESSOR: How does this help you match your read? How does this help you match your read, was the question. What is was your name?

AUDIENCE: Dan.

PROFESSOR: That's right, Dan. That's a great question. I'm so glad you asked it. First we'll get to Eric's question and then we'll get to yours. Because I know if I don't give you a good answer that you're going to be very mad, right? OK?

Let's talk about the question of how to do this without the other things. So we're going to create something called the last to first function that maps a character in the last row, column, I should say, to the first column. And there is the function right there. It's called LF.

You give it a row number. The rows are zero origined. And you give it a character and it tells you what the corresponding place is in the first column. And it has two components. The first is Occ, which tells you how many characters are smaller than that character lexically. So tells you where, for example, the a's start, the c's start, or the g's start.

So in this case, for example Occ of c is 4. That is the c's start at 0, 1,2, 3, the fourth row. OK? And then count tells you the rank of c minus 1. So it's going to essentially count how many times c occurs before the c at the row you're pointing at. In this case, the answer is 1 and you add 1 and that gets you to 5, which is this row.

So this c2 maps here to c2 as we already discussed. So this LF function is a way to map from the last row to the first row. And we need to have two components. So we need to know Occ, which is very trivial to compute. There are only five elements, one for each base and one for the end of line terminator, which is actually zero. So it will only have integers and count, which is going to tell us the rank in the BWT transform and we'll talk about how to do that presently.

OK. So did that answer your question, how to do this without the rest of the matrix? Eric?

AUDIENCE: Could you show us step by step on the blackboard how you would reconstruct it?

PROFESSOR: How do we reconstruct it?

AUDIENCE: Yeah.

PROFESSOR: You mean something like this? Is this what you're suggesting?

AUDIENCE: Somehow I get a feeling that the first column doesn't help us in understanding how the algorithm work only using the last column.

PROFESSOR: OK. Your comment, Eric, is that you feel like the first column doesn't help us understand how the algorithm works, only using the last column, right? OK.

AUDIENCE: [INAUDIBLE] going back to the first column of data.

PROFESSOR: OK. Well let's compute the LF function of the character and the row for each one of these things, OK? And that might help you, all right? Because that's the central part of being able to reverse this transform. So this is, to be more clear, I'll make it more explicit. This is LF of I and BWT of i, where i goes from 0 to 6. So what is that value for this one right here? Anybody know? Well it would be Occ of g, which is 6, right? Plus count of of 6 n g, which is going to be 0.

Or I can look right over here and see that in fact it's 6, right? Because this occurrence of g1 is right here. So this LF value is, it's 6 4 0 a1 is in 1, a2 is in 2, a3 is in 3, c2 is in 5. So this is the LF function, 6 4 0 1 2 3 5. And I don't need any of this to compute it. Because it simply is equal to, going back one slide, it's equal to Occ of c plus count. So it' going to be equal to where that particular character starts on the left hand side and its rank minus 1.

And so these are the values for LF. This is what I need to be able to take this string and recompute the original string. If I can compute this, I don't need any of that. And to compute this, I need two things. I need Occ and I need count. All right? Now, I can tell you're not quite completely satisfied yet. So maybe you can ask me another question and it would be very helpful to me.

AUDIENCE: How did you get G1's last and first functions score being 6?

PROFESSOR: OK. Let's take that apart. We want to know what LF of 6 and where was that G1? G1 is 1 and 0, right? Sorry. LF of 1 and g is equal to, right? Is that g and 1 or 0? Oop, sorry it's in 0. So this is what you like me to compute, right?

OK what's Occ of g? It's how many characters are less than g in the original string? I'll give you a clue. It's 1, 2, 3, 4, 5, 6.

AUDIENCE: [INAUDIBLE].

PROFESSOR: No, it's how many characters are less than g in the original string. How many things are going to distort underneath it? Where do the g's begin in the sorted version? The g's begin in row 6. OK? So OCC of g is 6. Is that-- are you getting hung up on that point?

AUDIENCE: Yes. How do you know that without ever referencing back to the first 5 columns?

PROFESSOR: Because when we build the index we remember.

AUDIENCE: Oh, OK.

PROFESSOR: So we have to, I haven't told you this, but I need to compute, I need to remember ways to compute Occ and count really quickly. Occ is represented by four values only. They're where the a's start, the c's start, the t's start, and the g's start. That's all I need to know, four integers. OK? Are you happy with that?

Occ, the a's start at 1. The c's start at 4 and the g's start at 6. OK? That's all I need to remember. But I precompute that. OK? Remember it. Are you happy with that no?

AUDIENCE: Yes.

PROFESSOR: OK. And then this business over here of count of zero and g. Right? Which is how many g's, what's the rank of this g in the right hand side? 1. Minus 1 is 0. So that's 0. That's how we computed it. OK?

These are great questions because I think they're foundational. Yeah.

AUDIENCE: Occ and count are both precomputed as you're building on this last column [INAUDIBLE].

PROFESSOR: They are precomputed and well, I have not told you, see, you're sort of ahead of things a little bit in that I'd hoped to suspend disbelief and that I could actually build these very efficiently. But yes, they're built at the same time as the index is built. OK? But yes.

OK and if I have Occ and count and I have the string, then I can get this. But I can still need to get to Dan's question because he has been very patient over there. He wants to know how I can use this to find sequence region of genome. And he's just being so good over there. I really appreciate that, Dan. Thanks so much for that. Are you happier? OK you're good. All right.

So and this is what we just did. The walk left algorithm actually inverts the BWT by using this function to walk left and using that Python code up there, you can actually reconstruct the original string you started with. So it's just very simple. And we went through it on the board. Are there any questions about it at all?

AUDIENCE: Yes. Can you actually do this any column? Like, why are you using the last column instead of like, why can't you just change it like the [INAUDIBLE], the equation and make it work for--

PROFESSOR: Because a very important thing is that this is actually a very important property right? Which is all the suffixes are sorted here. And if we didn't do that though, I couldn't answer Dan's question. And he'd be very upset. So please be respectful of his interest here. Now, the thing is, the beauty of this is, is that I have all these suffixes sorted and what you're about to see is the most amazing thing, which is that we're going to snap our fingers and, bang, we can map 200 million reads in no time at all.

You like that? You're laughing. Oh, oh, that's not a good That's not a good sign.

Let's press ahead fearlessly, OK? And talk about how we're going to use this to map read. So we're going to figure out how to use this index and this transform to rapidly aligned reads to a reference genome. And we're not talking about one read or 10 reads or 1 million reads. We're talking about hundreds of millions of reads. So it has be very efficient indeed, OK?

So here's the essential idea. There's the core algorithm on the slide. Which is that what we do is we take the original query that we have the read that we're trying to match and we're going to process it backwards from the end of the read forwards. And we begin by considering all possible suffixes from row zero to, in this case, it would be row seven. Which is the length of the entire transform.

And we iterate and in each iteration we consider suffixes that match the query. So in the first step, right here, let me see if I can get my point working, there we are. So in the first step here, we matching this c. OK? And we compute the LF of the top, which is this row and of the bottom, which is down here pointing off the end, and that takes us to the first d here and to this point.

Here are the two c's that could be possible matches to our query, which ends in a c. We then say, oh, the next character we have to match is an a. So we look here at the a we need to match, and starting from this row, which is row four, and this row, which is row six, we compute the LF of each one of these to figure out what rows in a precedes these c's.

And the way we compute the LF is that we use the character a to be able to figure out which rows have the a preceding the c. You can see, when we compute those LF functions, what we wind up with are these rows where we have a followed by c. So we're beginning to match the end of our read, as we go from right to left.

We then compute the same thing once again, considering the first a and ask what rows are going to allow us to put this a in front of the ac to form our entire read. And we compute the LF once again of these things. And you can see that here it takes us to this specific row aac. So that row represents a suffix that is matching our query exactly.

So we iterate this loop to be able to match a read against the index. And we're using the LF function to do it. And it's a really beautiful algorithm. And remember, we only have the transform. We don't have the rest of this matrix.

So before I press ahead and talk about other details, I think it's important to observe a couple of things that are a little bit counterintuitive about this. One counterintuitive aspect of it is, that when I'm over here for example, and for example when I'm computing the LF here, I'm computing the LF of row two with respect to a. But there's a dollar sign there. Right?

So I'm using this to the LF function, to tell me where a suffix would be that actually follows my constraint of having to have an a be the prefix of ac, where I am right now. This code is actually not fake code. It's the actual code that's in a matcher, for matching a read against the index.

Now let me just stop right here for a second and see if there any other questions. Dan is getting now his answer to his question, right? About how you actually use this for matching reads. You do this once for every read. And it is linear time. Right? It's the length of the read itself is all the time it takes to match in a huge genome.

So once we've built the index of the genome, in fact, most of the time when you're doing this sort of mapping, you don't build the index. You download the index off of a website. And so you don't have to pay for the time to build this index. You just download the index and you take your reads and the time to match all of your sequencing reads against a certain build of whatever genome you're using is simply linear in the number of bases you have. Questions? Yes. And say your name and the question, please.

AUDIENCE: How did [INAUDIBLE] come up with intuition [INAUDIBLE]? It seems like they just pulled it out of a hat.

PROFESSOR: You know, I asked Mike that the other day. I saw him in a meeting and he sort of surprised at how this has taken off. And he told me some other interesting facts about this, which you probably could deduce. Which is that if you only want to match reads that are four long, you only have to sort this matrix by the first four characters. But there are other little tricks you can play here.

Any other questions? Yes.

AUDIENCE: I'm Deborah. What is the FM index?

PROFESSOR: What is the FM index. Well, the guys who thought this up have the last initials of F and M, but that's not what it stands for, contrary to popular opinion. It stands for full text minute size. That's what they claim. So if you hear people talking about full text minute size indices, or FM indices, the Fm index is actually the part that was being asked over here, the Occ part and the LF part, how you actually compute those quickly.

That was what FNM contributed to this but, generically when we're talking about this style of indexing, it's called FM indexing or you might hear, I'm using a BWT. Some people will say that. But that's what FM stands for. Does that answer your question?

Excellent. OK. All right. Any-- these are all great questions. Yes.

AUDIENCE: [INAUDIBLE].

PROFESSOR: Oh, you don't know that a and c are there, except that remember, if you look at the way that this is working, is that you're not actually reconstructing strings, you're only trying to find them. Right? And so at the end, top and bottom are going to point to the row that contains the suffix where your original read was.

And now your next question is going to be, where is that in the genome? This doesn't do me any good. I mean, the number 1 doesn't help me out here, doesn't mean anything. Not good, right? So where is it in the genome is the next question.

So we'll get to that in a second. What happens if you give me a read that doesn't match anywhere in this index? Well if you give me a read that doesn't match anywhere in this index, what happens is the top and bottom become the same. So on top and bottom become the same, it's a failed look up. All right?

And that's because the suffix doesn't exist in the index. And once top and bottom become the same, they remain the same throughout that loop. Yes.

AUDIENCE: I'm Sally. My main question is that this doesn't provide any leeway for errors. You kind of have to be able to present all of your rates.

PROFESSOR: Sally, you're absolutely correct. I'm so glad you asked that question. Your observation is it does not provide any leeway for mismatches. And so unlike all the other algorithms we study, which had these very nice matrices and ability to assign weights to mismatches, this is only doing exact matching.

And so what you need help understanding is, how we can deal with mismatches in the presence of this. And I'll get to that in less than 10 minutes. And it won't be quite as elegant as what you saw from Professor Berg but it's what everybody does. So that's my only excuse for it OK? Yes.

AUDIENCE: What is the bottom set of arrows doing? What's its significance?

PROFESSOR: The significance of top and bottom, that's a great question. What the significance of top and bottom? Top and bottom bracket in that original matrix, the suffixes that are matching the original query. OK? And so between top and bottom minus 1 are all of the rows that have a suffix that match our original query. And if top equals bottom, there are no matching suffixes. But assuming there are matching suffixes, those are rows that contain a matching suffix. And as we progress along, top and bottom change as the suffixes change as we expand the suffix to contain more and more bases.

OK? OK. Any other questions? OK. So now back to the question over here, which is that OK, I know that we've matched, and I know that we have this hit. The question is, where is it in the Genome because the fact that it matched row one of my BWT matrix does me absolutely no good at all. Right?

Anybody have any ideas about how we could figure out where it is the genome? Given what I've told you so far, which is that you have the BWT, you have Occ, and you have count, and you can compute the LF function. Any ideas? Doesn't matter how slow it is.

OK well how could I figure out what came before aac in the genome? Yes.

AUDIENCE: So, for, like at the beginning, we rebuilt this whole string, starting at the end. You could rebuild the whole string starting at, rebuild the whole genome starting at the 1 and see--

PROFESSOR: You could rebuild the entire genome that prepends or occurs before aac, right?

AUDIENCE: Exactly.

PROFESSOR: Exactly. So that's what we can do. We could actually do our walk left algorithm. We can walk left from there, find out that we go two steps until we hit dollar sign, and therefore, the offset is two, where it occurs in the genome. So we can give a match position by walking left. Does everybody see that, that we could walk left to figure where it is? It's not fast, but it works. Yes.

AUDIENCE: I'm Ted.

PROFESSOR: Hi, Ted.

AUDIENCE: So now our function first has to take the read and it has to align it and the same, where the position is built into the end of the function, the speed of the function is now dependent on position as well. Is that right? Because the longer it takes,

PROFESSOR: I was being a little bit glib find if it matches or not this linear time. Now you're saying, hey, wait a minute. I want to know where it is in the genome. That's a big bonus, right? And so you'd like to know where that is? Yes. But I still can do that in linear time. And we'll show you how to do that in a second. This is not linear time.

This actually needs to walk back, order the length of genome for every single query. That's not good. All right? Well,

What we could do is, we could store what's called a suffix array with each row and say, where in the genome that position is. Where that row starts. And then maybe a simple look-up. That when you actually have a hit in row one, ah, OK, start your position two of the genome.

But then the problem with that is that it actually takes a lot of space. And we want to have compact indices. So the trick is what we do is, instead of storing that entire suffix array, we store every so many rows, like every 25 rows. And all we do is, we walk left until we hit a row that actually has the value, and then we add how many times we walked left, plus the value and we know where we are in the genome.

So we can sample the suffix array, and by sampling the suffix array, we cut down our storage hugely and it's still pretty efficient. Because what we can do is, we just walk left until we hit a sample suffix array location and then add the two numbers together. All right?

So that's how it's done. OK? So that's how we actually do the alignment and figure out where things are. The one thing I haven't told you about is how to compute count efficiently. Now remember what count does. Count is a function-- but this is putting it all together where we're matching this query, we do the steps. We get the match. Then we do walk left once and then we look at the suffix to figure out where we are, right?

The business about count is that what we need to do is to figure out the rank of a particular base in a position in the transform. And one way to do that is to go backwards to the whole transform, counting how many g;s occur before this one, and that's very expensive, to compute the rank of this particular g. Remember the rank is simply the number of g's that occur before this one in the BWT. Very simple metric.

So instead of doing that, what we can do is, we can build a data structure that every once in awhile, counts how many a's, c's, g's, and t's have occurred before now in the BWT. And so we're going to sample this with these checkpoints, and then when you want to compute count at any point, you can go to the nearest checkpoint, wherever that is, and make an adjustment by counting the number of characters between you and that checkpoint. Very straightforward. All Right

So this, coming back to question, it's Time, right? --asked you need to build this checkpointing mechanism at the same time you build the index, as well as the sampling of the suffix array. So a full index consists of the transform itself, which is the genome transformed into its BWT. And they literally take the entire genome and do this.

Typically they'll put dollar signs between the chromosomes. So they'll transform the whole thing. It takes a sampling of the suffix array we just saw and it takes the checkpointing of the LF function to make a constant time. And that's what is inside of an FM index. OK?

Now it's small, which is one of the nice things, compared to things like suffix tree, suffix arrays, or even other kinds of hash structures for looking for seeds, it really is not even twice the size of the genome. So it's a very compact index that is very, very efficient. And so it's a wonderful data structure for doing what we're doing, except we have not dealt with mismatches yet, right?

And so once again, I want to put a plug in for BWA which is really a marvelous aligner. And we'll talk about tomorrow in recitation if you want to know all of the details of what it actually takes to make this work in practice. Now, it finds exactness matches quickly, but it doesn't really have any allowances for mismatches. And the way that bow tie and other aligners deal with this, and they're all pretty consistent, is in the following way, which is that they do backtracking.

Here's the idea. You try and match something or match a read and you get to a particular point in the read, and you can't go any further. Top is equal to bottom. So you know that there's no suffix in the genome that matches your query. So what do you do?

Well, what you can do is you can try all of the different bases at that position besides the one you tried to see whether it matches or not. I can see the horror coming over people. Oh, no, not backtracking, not that. But sometimes it actually works.

And just to give you order of magnitude idea about how this works in practice, when reads don't match, they limit backtracking to about 125 times in these aligners. so they try pretty hard to actually match things. And yes, it is true that even with this backtracking, it's still a great approach. And sometimes the first thing you try doesn't work, and you have to backtrack, trying multiple bases at that location until you get one that matches. And then you can proceed. OK And you eventually wind up at the alignment you see in the lower right hand corner, where you're substituting a g for an a, an a for a g, excuse me, to make it go forward.

Do people understand the essential idea of this idea of backtracking? Does anybody have any comments or questions about it? Like ew, or ideas? Yes.

AUDIENCE: What about gaps?

PROFESSOR: What about gaps? BWA, I believe, processes gaps. But gaps are much, much less likely than missed bases. The other thing is that if you're doing a sequencing library, and you have a read that actually has a gap in it, it's probably the case you have another read that doesn't. For the same sequence. So it is less important to process gaps than it is to process differences.

The reason is that differences mean that it might be a difference of an allele. In other words, it might be that your base is different than the reference genome. Indels are also possible. And there are different strategies of dealing with those. That would be a great question for Hang tomorrow about gaps. Because he can tell you in practice what they do. And we'll get into a little bit of that at the end of today's lecture. Yes. Question?

AUDIENCE: I'm Levy.

PROFESSOR: Hi, Levy.

AUDIENCE: How do you make sure when you're backtracking that you end up with the best possible match? Do you just go down the first--

PROFESSOR: The questions is how do you guarantee you wind up with the best possible match? The short answer is that you don't. There's a longer answer, which we're about to get to, about how we try to approximate that. And what judgment you would use to get to what we would think is a practically good match. OK? But in terms of theoretically optimal, the answer is, it doesn't attempt to do that. That's a good question. Yes.

AUDIENCE: In practice, does this backtracking method at the same time as you're computing the matches or--

PROFESSOR: Yes. So what's happening is, you remember that loop where we're going around, where we were moving the top and bottom pointers. If you get to a point where they come together, then you would at that point, begin backtracking and try different bases. And if you look, I posted the BWA paper on the Stellar website. And if you look in one of the figures, the algorithm is there, And you'll actually see, if you can deconvolute what's going on, that inside the loop, it's actually doing exactly that. Yes. Question.

AUDIENCE: In practice, is the number of errors is small, would it make sense just to use [INAUDIBLE]?

PROFESSOR: We're going to get to that. I think the question was, if the number of errors is small, would it be good to actually use a different algorithm, drop into a different algorithm? So there are algorithms on FM index assisted Smith Waterman, for example. Where you get to the neighborhood by a fast technique and then you do a full search, using a more in depth principles methodology, right? And so there's some papers I have in the slides here that I referenced that do exactly that.

These are all great questions. OK. Yes.

AUDIENCE: If you're-- If you only decide to backtrack a certain number of times, like 100 times, then wouldn't like the alignment be biased towards the end of the short read?

PROFESSOR: I am so glad you asked this question. The question is, and what was your name again?

AUDIENCE: Kevin.

PROFESSOR: Kevin. Kevin asks, gee, if you're matching from the right to the left, and you're doing backtracking, isn't this going to be biased towards the right into the read, in some sense, right? Because if the right into the read doesn't match, then you're going to give up, right? In fact, what we know is the left end of the read is the better end of the read. Because sequences are done five prime to three prime and thus typically, the highest quality scores or the best quality scores are in the left hand side of the read.

So do you have any idea about how you would cope with that?

AUDIENCE: You could just reverse one of them. But you'd reverse--

PROFESSOR: Exactly what they do. They execute the entire genome and they reverse it and then they index that. And so, when they create what's called a mirror index, they just reverse the entire genome, and now you can match left to right, as opposed to right to left. Pretty cool, huh? Yeah.

So backtracking, just note that there are different alignments that can occur across different backtracking paths. And this is not optimal in any sense. And to your question about how you actually go about picking a backtracking strategy, assuming we're matching from right to left again for a moment, what you can do is, if you hit a mismatch, you backtrack to the lowest quality based position, according to PHRED scores.

We talked about PHRED scores earlier, which are shown here on the slide. And you backtrack there and you try a different base and you move forward from there. So you're assuming that the read, which is the query, which is associated quality scores, is most suspect where the quality score is the lowest. So you backtrack to the right to the leftmost lowest quality score.

Now it's a very simple approach. Right? And we talked a little bit about the idea that you don't necessarily want to match from the right side and thus, typically the parameters to algorithms like this include, how many mismatches are allowed in the first L bases on the left end, the sum of the mismatch qualities you're going to tolerate, and so forth. And you'll find that these align yourself with a lot of switches that you can set. And you can consult with your colleagues about how to set switches, because it depends upon the particular type of data you're aligning, the length of the reads and so forth.

But suffice it to say, when you're doing this, typically we create these mirror indices that actually reverse the entire genome and then index it. So we can either match either right to left or left to right. And so for example, if you have a mirror index, and you only tolerate up to two errors, then you know that either, you're going to get the first half right in one index or the mirror index. And so you can use both indices in parallel, the forward and the reverse index of the genome, and then get pretty far into the read before you have to start backtracking.

There are all these sorts of techniques, shall we say, to actually overcome some the limitations of backtracking. Any questions about backtracking at all? Yes.

AUDIENCE: Is it trivial knowing the BWT originally to find the mirror BWT? Like for example,

PROFESSOR: No, it's not trivial.

AUDIENCE: So it's not like a simple matrix transforming [INAUDIBLE].

PROFESSOR: No. Not to my knowledge. I think you start with the original genome, you reverse it and then you compute the BWT with that. Right? That's pretty easy to do. And Hang was explaining to me today how you compute his new ways of compute BWT, which don't actually involve sorting the entire thing. There are insertion ways of computing the BWT that are very [INAUDIBLE], and you could ask him this question tomorrow if you care to come.

All right just to give you an idea on how complex things are, to build an index like this, takes, for the entire human genome, we're talking five hours of compute time to compute an index to give you an order of magnitude time for how to compute the BWT. the LF checkpoints, and the suffix array sampling. Something like that.

So it's really not too bad to compute the index of the entire genome. And to do searches, you know, we're talking about, like on a four processor machine, we're talking about maybe upwards of 100 million rads per hour to map. So if you have 200 million reads and you want to map them to a genome, or align them as it's sometimes called, it's going to take you a couple hours to do it. So this is sort of the order of magnitude of the time required to do these sorts of functions.

And there are a couple fine points I wanted to end with today. The first is we haven't talked at all about paired, erred, and read alignment. In paired read alignment, you get, for each molecule, you get two reads. One starting at the five prime end on one side, , and one starting from the five prime end on the other side. So typical read links might be 100 base pairs on the left and 100 base pairs on the right.

What is called the insert size is the total size of the molecule from five prime end to five prime end to read. And the stuff in the middle is not observed. We actually don't know what it is. And we also don't know how long it is. Now when these libraries are prepared, size selection is done, so we get a rough idea of what it should be. We can actually compute by looking at where things align on the genome, what it actually is. But we don't know absolutely.

If we were able to strictly control the length of the unobserved part, which is almost impossible to do, then we would get molecular rulers. And we would know exactly down to the base, whether or not there were indels between the left read and the right read when we did the alignment. We actually don't have that today.

The sequencing instrument actually identifies the read pairs in its output. That's the only way to do this. So when you get an output file, like a fast Q file, from a sequencing instrument, it will tell you, for a given molecule, here's the left read and here's the right read. Although left and right are really sort of misnomers because there really is no left and right, right? This is one end and then this is the other end.

Typical ways of processing these paired reads, first you align left and right reads. And they could really only be oriented with respect to a genome sequence where you say that one has a lower coordinate than the other one when you're actually doing the alignment. And if one read fails to align uniquely, then what you can do is, you know what neighborhood you're in because you know, roughly speaking, what the insert size is, so you can do Smith Waterman to actually try and locate the other read in that neighborhood. Or you can tolerate multiply mapped reads.

One thing that I did not mention to you explicitly, is that when you match the entire query, and top and bottom are more than one away from each other, that means you've got many places in the genome that things map. And thus you may report all of those locations or I might report the first one. So that's one bit of insight into how to do a map paired reads.

And these are becoming very important because as sequencing costs go down, people are doing more and more paired and sequencing because they give you much more information about the original library you created and for certain protocols can allow you to localize events in the genome far more accurately.

Final piece of advice on considerations for read alignment. We talked about the idea that some reads will map or align uniquely to the genome and some will multimap. You know that the genome is roughly 50% repeat sequence. And thus it's likely that if you have a particular read molecule, there's a reasonable chance that it will map to multiple locations. Is there a question here? No. OK.

You have to figure out what your desired mismatch tolerance is when you're doing alignment and set the parameters to your aligner carefully, after reading the documentation thoroughly, because as you could tell, there's no beautiful matrix formulation like there is with a well established basis in the literature, rather it's more ad hoc. And you need to figure out what the desired processing is for paired reads.

So what we've talked about today is we started off talking about library complexity and the idea that when we get a bunch of reads from a sequencer, we can use that collection of reads to estimate the complexity of our original library and whether or not something went wrong in the biological processing that we were doing. Assuming it's a good set of reads, we need to figure out where they align to the genome.

So we talked about this idea of creating a full text minute size index, which involves a Burrows-Wheeler transform. And we saw how we can compute that and throw away almost everything else except for the BWT itself, the suffix array checkpoints, and the FM index checkpoints to be able to reconstruct this at a relatively modest increase in size over the genome itself and do this very, very rapid matching, albeit with more problematic matching of mismatches. And then we turned to the question of how to deal with those mismatches with backtracking and some fine points on paired end alignment.

So that is the end of today's lecture. On Tuesday of next week, we'll talk about how to actually construct a reference genome, which is a really neat thing to be able to do, take a whole bunch of reads, put the puzzle back together again. I would encourage you to make sure you understand how this indexing strategy works. Look at the slides. Feel free to ask any of us.

Thanks so much for your attention. Welcome back. Have a great weekend. We'll see you next Tuesday.