Video 2-9: Motifs finding heuristics

0:02 The rest of what I want to tell you maybe also is a a short story. Um I want to present some some maybe a generalizable recipes for doing these uh these heristics. So um solving a problem that's pro that's that's provably hard where no obvious thing works and uh we decide to give up the optimality guarantee. So we're willing to take some solution and uh essentially hope and pray that it's a good one or check in in with other means and um give up provable optimality. Uh that happens more often than not in in some binformatics problems.

0:47 Now yeah uh terms of uristic can mean many things to many people. In the context of algorithms a uristic is basically nothing. I mean it's uh it's an algorithm that produces a valid solution to your problem. It's not going to just uh give you some some illegal output. But the guarantee can there's nothing known about the guarantee of the the quality of the output. Uh it's to be contrasted to approximation where an approximation algorithm has to have a specific guarantee. For example, C approximation means at most a factor C worse than the optimum. Uh but uristic uh cynically you could say is uh some person has one input where it worked and that's as much as we promise. Um it yeah if if you had a bit more time maybe it would be fun to brainstorm what you would do on these um on this motive problem. You're given this collection of strings and how do you even get started with doing something that's at least not taking forever even if it might not find the best possible indices for this for the problem the best consensus string.

2:13 Uh a simple thing you could try um is to kind of go one string at a time and try to build uh um the most sensible continuation of selections of strings. So for the first string I'm going to go over all possible start positions. So I keep that brute force bit. Now, if I fix some part in the first string, should have done that in red, I guess.

2:55 Now I'm looking at the second string. How can I find something that is a sensible continuation of the red string? Um, for that, let me borrow the picture from earlier again.

3:16 One option or one one sensible option would be to just say uh what we want is u a region of that substring that doesn't have many um mismatches with the first first row So you could say um we're in this situation now that we've selected this string and now we're looking at the second we're trying to fill the second row but we have a much longer string and we have to decide which one to take. Uh things are easy if we already have as in this picture highlighted which characters are eventually going to be to the the contribution to the consensus string. And one option would be to say I'm going to compute the consensus of what I have so far and then find the best match for that consensus string. But it turns out that's a bit um that's a bit more blind than we can be because it ignores things where um there's a close tie between two characters and then a third character is a mismatch to both whereas one of the two characters might currently be in the minority but might soon turn into a majority and things like that.

4:47 So um what we'll do instead is look at uh the the profiles as it's as is uh explained as it's called here. If you have um suppose this is our our entire picture. So we've already selected a certain number of indices. We're further into this iteration not just uh at position two but we're now adding the next row from the next string. Uh how do we find that row? We compute the profile which is we sum up in each column separately how often each of the four characters occurs. That's what the profile or the well they call it count here in profile and is normalized to sum to one but of course the same information. What that gives you is uh a way of saying yeah there's uh a column could be completely messed up in the sense that all characters occur equally often or it could still be not a clear majority but it's only messed up between two characters or it could be like like this column predominantly a and you can quantify that that's what's meant by profile. So you take your current selection from the starting indices you have. You take the current selection of K mirror of of length K substrings. In each column you separately count how often the letters occur and you turn that into probabilities.

6:17 Uh I'll comment on the zeros in a second. So that's what's meant by profile here. Okay, that's this one. uh so there's a different profile for each position and what the profile is is a a probability for each letter then I can take such a profile as my stochastic model to uh evaluate a new string where I could generate random strings from it where each position is independently drawn from the distribution given by that that column's probabilities and so I can uh I can write this this conditional probability. Um maybe we do a tiny example just to make clear the definition is is uh is understood. If the profile looks like um so I have a cg t then I have the columns 0 1 2 I'm just writing the profile as an inline table um whatever encoding B really uses suppose uh here it's I don't know A has five C has 2.1 1 here it's 2.5 2.1 I'm only going to use the same numbers I don't mess up this the sum and the conditional probability of producing a a under this profile would be the entry in the first column for a.5 times the probability for a in the second column and then the c in the third column that would be.5 that's the the conditional probability of uh of generating this thing among all possible strings of that length. Okay.

8:36 And so here um what I'm doing for the next string that I pick my my window from, I iterate over all the possible start positions, take the length k substring and compute how likely is that to arise if I'm assuming the current profile from the previous rows as the truth. And then I just pick the one that's most likely. That's what what this does. It's a And then um I'm doing this over all the rows. So I keep growing my selection of substrings. Each time I'm picking the one that assuming the prefix I've computed so far is the truth is most likely uh among the options that that new string gives me. And if I've done this for all I've I've selected a certain collection of substrings. Now I can check if their score is better than what I previously had. And here this this um this just computes the up the maximum or the minimum over all the start positions in the first string.

9:54 Right? For each position in the first string, the rest is completely deterministic. And that also means the whole thing is deterministic. There's no randomness involved so far. Even though it uses these probabilities, it's just using them to make a selection. But I'm picking a deterministic maximum.

10:23 it's deterministic but uh it's very sensitive to the order in which I look at things and the order in which I list the genes. That's very arbitrary. And so uh it might work or it might not work and it depends on which is my first genome mostly or it also depends on the order of the others but the first one I'm brute forcing over the other I'm just doing greedy stuff so it's a big difference and that seems for the example um of finding this this gene in plants it seems to not work at all. There's another subtle problem. Um, which brings me back to this picture. If you do what I described so far, it certainly doesn't work well. Um, but the reason is stupid and fixable in the profile as we defined it. There's a ton of zeros.

11:20 If you look at at this or I mean the counts is the same. Uh, there's a lot of zero entries here. And because we're taking products, if you have any column that contributes a zero to your position in the string, your likelihood is zero. So especially initially, if you think about the first step, the first is you pick some substring in the first row, then you look at the second row, the profile is going to be 100% deterministic, right? There's only one character each. So the profile will say first column has to be an A, second has to be a Z, everything else is probability zero. So it's only going to pick uh either there's an exact batch in the second row, then you pick that one, or it's going to pick equally likely any of the others. That's not that useful.

12:18 So in a way, we're too strict. It's it's going in the right direction. This maximum likelihood idea seems uh sensible, but we're way too strict. And the problem is that something times 0 is exactly zero no matter what the other thing was. So we can fix that um in in short by starting to count at one instead of at zero. [snorts] There's many ways to justify this. Starting at one is not magic. You could start at some other number. All that it does is essentially you're mixing in a prior distribution that is uniform and then you're doing Beijian inference on when you see characters it moves away from the uniform but in a smooth way. If you don't start with one but with some other parameter your uni your prior is not uniform but some uh directly distribution with a different parameter. Uh I mean I I'll leave that fun to the statisticians. Uh for for us if you're collecting enough columns it doesn't matter all that much what initial value you put as as long as it's not identically zero.

13:54 So this is an important fix. If anything of this sort with using profile probabilities, profile driven probabilities, you have to have to have to start with something other than zero when counting.

14:11 Uh this is this is sometimes called a plus rule. I think it's somewhat folklore and goes by many names, but keep that in mind whenever things seem not to work. uh in in some kind of a statistical inference thing where where you expect non-exact answers maybe start counting from one is is already fixing it. Okay, there's this other problem that really this is so sensitive to the order and we have no good reason to prefer one order over the other that it's not clear how to really fix this one. Uh there's another approach that is very general. We start with something random and then we try to locally improve it. This one was a bit ad hoc. It kind of it started also with something quote random by picking any starting index and then we try to brute force over that first choice. But then the rest was kind of uh greedy incremental. Maybe we can do something else instead. We start with a random choice among the collections uh among the possible starting indices. And then for that collection we try to modify it to make it slightly better.

15:27 And I don't really know how know how to better put this. Um so let me just first give some some pictures and the code maybe is self-explanatory if you read it but I don't want you to first read first listen. So if you take some initially random selection of these of these things that's going to give you some consensus string I called it t bar earlier let's use that that's short now I can take that consensus string and go back into the strings and select the most likely red thing or that consensus string.

16:14 So from the starting indices I can go to some T bar and from those I can go to some other starting indices and then from those I can go back to a new T bar that might be different and I can keep going and in each step I'm only going to make the score exactly better or it stays the same. It can never go worse by definition. I'm I'm choosing the update from S to T bar so that it exactly is uh the best choice for T bar and that if I pick the optimal T bar that gives me exactly the definition of this DH and if I go from T bar and find a better placement of these red things that also is only going to reduce the cost.

17:08 And so in code it looks it looks more complicated than um it really is. Uh you start with a random initial selection here you compute the profile um then you pick for that profile instead of using the consensus string uh is more robust to use the profile again. you could use the consensus string, but by the same reasoning, consensus is this hard maximum. It's a bit more robust to use the profile. Um, then you pick the most likely starting position given the profile. If that's better, then you use that instead. And if that hasn't improved the current solution, you're actually stuck. It's not going to work if you try again. And then you just stop. That's why it's a a hill climbing method. It finds a local optimum. You have no guarantee if that local optimum is any good compared to the global optimum. Um I guess it should it should climb down because we're minimizing but uh ignore that kind of hill climbing is always formulated as maximization.

18:17 You're rolling down into the steep the local valley. That's another option. uh once you selected this random starting position, it's deterministic and it's always going to reach some local optimum. Just no idea how good that is. Um in general, that method if you just run it once is not performing well. What you can do is run it many times and give a a certain a certain balance between how long you keep trying to make improvements. Maybe you stop it early because if you're in a very bad space, uh local improvements are not going to cut it and you just stop it and run it again with a new random starting point. Um, so you would probably want uh because of that you want to repeat a few times.

19:40 uh but unfortunately we don't have an analysis of how to choose these parameters how often to try for how long to try um unlike maybe things like shurning's threes algorithm where you could precisely analyze how long you have to run it so that it's in expectation finding the right answer. So there's in general nothing known as far as I know.

20:06 Okay, we're out of time with the planned five minute stop early. There's there's one last method. Uh maybe we'll do that next week then, unless you want to see the the final answer from the book.

20:22 It's a twist of this one. Basically, this uh this arc max can be weakened a bit. There you go. uh you also start with a random initial um choice. Uh here it's formulated to repeat a specific number of times because this is not necessarily reaching an optimum. It keeps trying. Um and this time we're just replacing one of the red collections. So we have our current selection of uh substrings. Then we're randomly selecting one of them, this one. And now we're going to get rid of of the old one and choose a new one for that, but just for this single one. And that turns out to be a bit more robust. Oops. Stay where you are.

21:19 So, we're going to update that to a new one. Maybe over here. And how do you do that? Uh, you compute the profile for everything without that row. And then you could have done the same as before. Take the most likely new position. Um a twist here with this that's then called gibs sampling. Um you you compute for every possible position in this in this string where we're selecting a new thing from. We're computing for each possible position the probability that it arises under the profile. And then we sample with probability proportional to that. So previously we would have said a hard max find the most likely position take that now we're saying let's do a random choice and the probability for picking each position is the likelihood under this under the current profile. So that means it can also go back to worse pos worse solutions with a certain probability and that's why it's not always terminating but then it also has the chance to escape local optima and yeah the rest is I I think um same as before this method seems to be working best for this kind of problem in the applications where people use the for modifying is data dependent but it seems uh This choice is robuster is more robust than taking the maximum even if it's a bit slower.