Video 4-3: de Bruijn graph assembly, Hierholzer's algorithm

0:03 maybe somewhat slightly surprisingly at least uh even though this overlap graph leads to an NP hard problem Hamiltonian path that comes out of our uh sequencing a se sequence assembly problem there's another way to look at the same problem that does allow a poly time solution and actually a very efficient one uh let me tell a little bit of the backside story and this is I think 194 40s is the right time frame for this. Um I'm going to butcher the uh the Dutch name Nicolas Go De Bruyne. Um the last name is called is is pronounced De Bruyne I learned but uh I'm not sure about the first name. Uh he he he puzzled he wondered over pondered over the following he pondered the following puzzle. Suppose you want a single string that contains all possible length k strings over zeros and ones as a substring. How long would that string have to be?

1:05 Now why he wondered about this is a different question. [laughter] Fair enough. It's a math puzzle. And I mean as long certainly there is such a string. You can just concatenate all the 2 to the k strings. So upper bound very clearly 2 to the k * k right all one after the other. That's hardly going to be the shortest because they're all using the same characters over and over. You can reuse characters massively, but that's an upper bound. So it's certainly possible to do this. The question is what's the shortest and how can you find one, right? Um, as any good computer scientist, if you see a new good puzzle, you try out small values. So, here's the smallest sensible value, I guess. Uh, for K equals 2, there's four different strings that we want to see. And now let's do some some uh lower upper bound business. Uh, let's do some lower bound.

2:09 This uh two universal string, that's how he called those. A string is universal if it contains all the possible cameras. Uh it needs to it needs to have at least five letters. Why? Every possible of the four needs to have a different starting position because they're all different. And so I need four positions where things can start. But then from the last one, I I'll have one character stick out, right? So, I need at least five letters to make a two universal string over the binary alphabet. Is that clear?

2:45 Um, turns out you can also do this. Here's one. [laughter] So, you have 0 0, you have 1 one, you have 0 1, you have 1 zero. All four occur. And we've used five characters. So, clearly this is an optimal solution because I've just proved everything has to use at least five. curious h but k equals 2 is is tiny right so uh maybe you have to try another one uh here's for k equals 3 now you can do the same kind of calculation you'll have eight different characters they need different starting positions and the last one will have two characters stick out so that's 10 and this one has 10 characters and you can check that every possible possible three bit string occurs as a substring. So again we could solve it optimally like a smallest possible length. Now maybe uh you feel there should be a pattern right? Okay, that's what De Bruyne wondered about and he eventually proved that you can always find so there always is a shortest K universal string that is as short as can possibly be and you can find one efficiently and this is true not just for the binary alphabet but for any alphabet and for strings or for ks of arbitrary size. So just doing examples wouldn't get you there. But the proof technique is very elegant and it leads to a solution to the sequence assembly problem. So what did he what did he do?

4:24 Um yeah I said this uh he constructed a different graph and that's now called the de bruin graph. It starts also with our uh reads. So our kert composition. Now this time we'll use uh again vertices that correspond to strings but there's two differences. These are only k minus one letters each and I'm uh I'm going to compare uh kind of fold together uh duplicate elements. So if there's several reads uh several um kamers in the input that have the same first k minus one letters and that's especially true if there's duplicates in the input. Uh I'm I'm uh treating all of these prefixes only as a single vertex.

5:20 The clever part are the edges. And I'll have an edge for every for every kr in the input. And the edge connects the vertices that correspond to the first and second k minus one characters. So here's say we pick some example. Um this edge corresponding to CCA that's one of the kmers in the input. Three mers in the input. It starts with CC and ends with CA. So, it's an edge that goes from CC to CA.

5:52 Okay, let me maybe copy this here just to have that explicit. So, all the crucial information really is in the edges this time and the vertices are just there to connect things up. Okay. Now, uh I should have should have asked you this, but it's already on the slide. Uh in analogy to what we've previously had there, the vertices were in correspondence with the kamers, and I needed to find a path that uses every vertex. Now, the kurs are on the edges. So, I want a path through this graph that uses every edge once, exactly once. And that's called an oiler path. I'll I'll introduce this a little bit as a recap if if you've forgotten. Um the the main thing is oiler paths unlike Hamiltonian paths it sounds like a minuscule difference. Visit every vertex, visit every edge. There's a massive difference in terms of complexity. You can compute these very efficiently. And so we've transformed our problem into a different graph problem. We're still solving the very same optim the same uh string problem. But now suddenly there's no NP heart problem involved.

7:14 Of course as a digression this this kind of means maybe that also the special inputs we get in the Hamiltonian path problem that result from you from the overlap graphs may not be uh a hard subset of all possible Hamiltonian path inputs. Maybe that's also a way to make progress. Uh but that would at least need studying carefully what these graphs look like. Here we don't even have to whatever the de bruin graph looks like. I can always find oilarian paths.

7:46 All right, here's um an example again. Um let me highlight one or two things. As before um in this in this case where uh the start and end are unique to mayors. They don't occur anywhere else. You have vertices that are uniquely the start and end point of any h of any uh oiler path.

8:09 There's three occurrences of ATG, the blue ones. And so there's three edges for ATG. Of course, there are parallel edges. So this is a a multigraph where I might have several edges connecting the same pair of vertices. Of course, these have to be parallel edges because they're the same strings. They'll go from the same source to the same target if I take substrings. So that shouldn't be surprising. Um what else? Uh we still have this infamous self loop from uh this uh triple G, but this time it's crucial that we include it. So multigraph here also means allowing self loops. And it's important because we want to take all edges right now. Of course, the self loops are kind of for free. Whenever you visit this vertex, you just do the excursion and then you're done with this edge. Okay?

9:05 Now we can use these de Bruyne graphs also to solve the original puzzle that De Bruyne started with. You can build uh this graph from the set of all possible kares say all possible binary kars or all uh over a general alphabet then you'll have uh sigma to the k minus one vertices because all possible k minus one substrings occur at some point and for every possible sigma to the kr we have one edge.

9:45 Now this needs a little bit of inspection but it's not too hard to show just by symmetry basically that every vertex has the same in andout degree namely sigma and maybe just from like a small example suppose k is three and sigma is two so we're back to the binary alphabet and this might be one vertex what are the in and out uh edges of this of this vertex The out edges are all the three mares that I can form that start with 01. And that's these two 0 1 0 and 0 1 1, right? And they're going to lead to some vertex that's dictated by whatever is here, but I don't even care so much. Uh same for the incoming ones. They have to end with 01. Then I can have a zero or a one in front and that's all. So the crucial clever insight here is the outderee and the in degrees are equal and that will uh again I'll I'll recap this in a sec that will mean that the graph is oilarian. So there is an oilarian cycle and so if you follow that cycle you visit every s every of the sigma to the kmer exactly once but you've also just used that many steps. So the total length of this is going to be exactly sigma to the k characters plus k minus one and that matches our lower bounds. If you insert sigma 2 and k equals 2 uh then we're we're back to five. That's the example we've seen before. And 10 for k equals 3, right?

11:39 So in a way um the bronze puzzle is solved and we didn't even need to do much uh in this in this case all the observation came from symmetry.

11:49 The reason why this graph is balanced same inderee as outderee in every vertex was just out of symmetry of the construction. In the general case when we start with some k composition of a of a text it won't be the graph won't have all possible kmers k minus one mer as vertices and it may have uh more interesting in in andout degrees. So let's let's briefly revisit um what is known about these oiler oiler walks or oiler cycles. I'll use these mostly interchangeably. I should probably train myself to say oiler walk whenever I talk about these. A walk is just anything that follows edges, but it might reuse vertices as many times as we want. Um, now a walk in general can use edges several times as well. Here I really want um an edge simple walk. Every edge can occur at most once. And then it's an oiler walk if every edge occurs exactly once. Okay.

12:54 Uh the picture here is the old city of Kernixburg now uh what's it called uh Khalinrat um which used to have these bridges and that's where Oiler was first asked is there a way to walk across all of the bridges exactly once and come back where you started and well the story has it that's how Oiler kickstarted graph theory he represented the problem as a little graph and then he noticed nope In this graph, there is no oiler cycle because of course he didn't call it like that back then. Well, it was an interesting character apparently but not not that uh uh pompous. Uh you can easily prove this theorem. There's an oiler walk in your graph if and only if the number of vertices with odd degree is either zero or two.

13:51 If it's zero, then there's a cycle. If it's two, then there's a path that uses every edge. And if neither of the two applies, then there is no such thing. And you can prove both implications. Um, I'll I'll leave that leave that out for now because this is on an undirected graph. The original problem here, the uh bridges were not one-way streets. You could walk them in either direction and that would have would have been fine. And of course, you know, Oiler walks from kindergarten. uh and I'll further comment on this. How can you do this systematically? Here's an algorithm that um I always show in efficient algorithms as well. Uh there as an application of depth first search that's really what you can do. Um maybe let me just illustrate it briefly on this um on this what in Germany is known as house of Nicolaus.

14:56 Uh you can you can count the degrees and there's two vertices with odd degrees. So we'll have to start and end in those two. And the goal is to draw this this drawing in in one go without doing a jump with your pen. That's exactly an oiler walk. You can draw all the edges without lifting your pen. So we might start here.

15:17 And basically we're doing um a traversal of the graph. And whenever we get stuck we uh label the edges backwards. Okay. What what I mean by that? Um let's see that we get stuck. Let me walk here. So I'm just following edges blindly. Um, and so at this point I'm stuck. So now I'm going to count the total number of edges here is 1 2 3 4 5 6 7 8. Is that right? Um, so this edge that I last took where I got stuck, I'm going to label this with an eight and then backtrack over it as in a depth first search. Now at this node, uh, I could have taken another edge. So I'm doing that. Then I'm taking these and now I'm stuck here because I've used up all the edges. Now I go back over this one and use the next next smallest free number. Back over this one here. I came from here. You have to imagine doing a depth first search where you keep these on a stack.

16:25 So I'm remembering them now because it's it's a tiny graph. Uh if you properly do it, you'll have a stack of them. So here we came from down and then from here and two and one. And if I'm done with this, the numbers on the edges will form an order cycle. 1 2 3 4 5 6 7 8. Okay.

16:51 And it's not too hard to prove that this always works. Uh I'm not going to do it in in more detail here. The code is very detailed. Um the the key thing is it's essentially a depth first search. The difference is we're uh not stopping visiting a vertex once we've seen it for the first time. It's an edgebased depth for search. So I I will revisit vertices as long as there's edges I haven't seen. So um I have a a visited array which here is two-dimensional. It has an entry for every possible edge. This could be a map from edges to true false. Um it's just to write it write it more concisely.

17:39 Otherwise it's a standard depth for search. You have a stack of vertices where you came from. In each point we start with putting some start vertex on the stack. Here we have to make sure we use one of those odd degree vertices if there is one otherwise any vertex can be used as a starting vertex.

17:57 Um I'm remembering for every every vertex in the graph what is the next edge to uh to look at when I'm visiting that vertex again. So I'm remembering an an iterator in its adjacency list and uh initially that's just pointing to the first edge of every vertex. Then I'm calling the recursive uh depth first search. It's with an explicit stack. Sorry, nothing nothing recursive. Um it's the iterative version. Uh I have my start vertex that initially goes onto the stack. And then as long as there are things on the stack, I'm not done yet. I get the topmost vertex and look at the next edge from that vertex that is not yet processed.

18:44 Uh if there's if there's none left, that's the backtracking case. If there is one left, then um I'm following that edge. So uh there's again the case um sorry that's the case. I follow the edge. If it's if it has a next edge, then I'm following that edge. I mark the edge as visited and because it's an undirected graph, I mark it visited also in the back backwards direction.

19:12 And then I'm uh pushing the vertex I arrived at. This is the backtracking case. If a vertex I reach doesn't have any more outgoing edges, I have to go back from where I came. So I'm popping the last vertex from the stack. That's where I came from. And I'm writing in this oiler array uh at the last position j starts at m I'm putting in the uh uh the edge that I've uh used here that I've used to get here. Okay, that's the same algorithm I've showed in the example u just uh explicitly written with a stack as you would for a dep for search.

20:00 Um all of this was for undirected graphs. That's how you usually introduce Oilerian cycles. Uh for our de Bruyne graph, we need this in a directed graph which maybe sounds like we have to do all the work again, but there's only very minor changes that needed. First of all, Oilers's theorem um if you if you think about the proof that I skipped, the directions are uh both simple and if you take them separate is is very easy. If we have the second condition, then we've just listed an algorithm that computes an oiler walk. So certainly one exists. So [snorts] that direction is given by the algorithm. The other direction is easy. If it has an order walk, you can just observe that apart from maybe the end start and end point, every vertex in the middle is entered and left the same number of times. So it has to have an even degree in this undirected case. Right?

21:03 In the directed case, something similar is still true. If there is an oiler walk that is directed, well, every nonstar end vertex still has to be entered and left the same number of times. So, its total degree will have to be even. But actually, because the edges have directions, it has to have the same in-n-out degree. So, it's balanced.

21:28 The connection here is the only the only reason why this is more complicated is uh you can't concisely say the number of odd degree vertices is zero or two because for the start and end vertex you either have one more incoming or one more outgoing edge you need one of the two types. [snorts] So you can lift all the theory to directed graphs and the same kind of proof works and the same kind of algorithm works to compute the oiler walk if the right hand side is is fulfilled with um minor adjustments. Um initially when we started here we looked at the degrees of vertices and then picked one that has odd degree if there was one. Also if uh we're not in the case that the graph is Oilarian, we just stop.

22:22 Uh we can do the same here using this this condition. The start vertex has to be the one that has more more outgoing than incoming edges. Also obvious. So the only non-trivial change is at this point where we mark edges true where we mark that an edge is visited because it's directed. They're not the same. So we don't do this. This is left out in the directed case. All the rest of the algorithm is the same. That's the beauty of going once through the pain of writing down a dep search version. Then uh the directed case just pops out as unchanged.

23:07 Okay. Uh I went through this fairly quickly because I think everyone will have seen that at some point in their life but it might be a while ago. So uh also the directed version is often treated a bit stepmotherly like a [laughter] uh for some reason this is not receiving the same amount of love probably because of the historic stories.

23:29 Uh also it's it's fairly easy. Now let's do um a little bit of an uh a detour on the analysis of the algorithm. Uh how efficient is the algorithm? It's a depth for search even though it's it's edge centric. It essentially visits everything in the graph once but it visits every edge at least at at most once. But the running time otherwise is actually the same as a standard depth first search because that also visits every edge. It just doesn't traverse it and moves on. It just uh looks over the edge and says ah the target I've seen before. I'm not going to go there again. But you're still going to have running time proportional to the number of edges as well. Uh the same same is true here. You're essentially treat every edge once. Fair enough. maybe twice because you also uh mark it with an index. But a constant factor aside, the running time is the same. Now uh what's more expensive in this implementation is the space that we use at least in a graph where there's many more edges than vertices.

24:44 the um data structure for the stack here, this frontier of vertices that we're still active actively exploring, this may contain the same vertex many times. Um I don't know if this is uh if this is obvious from the quick pace I went over this, it is certainly obvious from the example we've done that we revisit vertices. So the condition here is I may have a vertex any number of times on the stack as long as the stack forms um an edge simple walk. So I'm not going to lift have a vertex more more copies of the vertex on the stack than it has edges. So in total the size of the stack is bounded by the number of edges but it may be much more than the number of vertices if there's more edges than vertices. Okay. So worst worst case you're walking a very long path and essentially using all the edges then all of the copies of the vertices are on the stack simultaneously.

25:54 And so as a as a fun little detour of some work we've done in our in my group um we asked ourselves is this is this needed? We're doing a depth for search even though a little modified one. Do we really need m space, mwords of space? Uh or can we get get rid of this dependence on the number of edges and go back to what standard traversals have m plus n time but n space because every vertex is at most once on the stack in the standard DFS.

26:26 And the answer is yes. That's a paper coming up in January on the symposium on simplicity in algorithms. So it's not just supposed to be um a new result but also a very simple algorithm. So I thought we might discuss it here uh briefly. Um and the key observation is that really for oiler walks you don't actually need the full power of a depth first search. If you think about um if you think about it this way again my uh tiny example from before there's another vertex. Um I can basically walk whatever when I'm stuck I'm just going to revisit some place where I can I can branch off and insert that into the current tour. So I could have said what if we do this one here I was stuck the red to is is done. Uh but I could have branched off here with a blue tour which is a cycle and it necessarily let me draw it a little little less uh confusing.

27:42 Um it has to be a cycle because after I've taken out the red one is a path. So it's taken out all the odd degree vertices. What's left is only even degrees. So any oiler paths, oiler walks, walks, walks, walks in the remaining graph will have to be cycles because there's only even degrees. And so I can take the red tour up until this point, then take the blue detour and then continue.

28:12 And that always works. Now, if I have several of these blue cycles, any combination, any ordering of inserting these other cycles into the red path will be a valid oiler walk. It doesn't matter. I don't have to remember in what order I found them. Now, a depth for search very carefully keeps track of the precise order in which I found all those edges. And the order walk just doesn't care. So, let's not do that. Let's not store the exact order in which we found edges. Um, but you still have to make sure you don't get stuck in the final backtrack. So, you need to store the first edge you reach a vertex with separately, but then all the others you can just say, "Yeah, I've seen those." The order doesn't matter. Uh, let me show on a on a short example how this works. And uh I well I couldn't convince my co-author my co-authors to convert this to something that I found more more intuitive. So we'll we'll do the whole algorithm backwards using all the edges backwards with the benefit that you can produce the oiler walk in the right order. Whereas here remember um here's algorithm here a while ago. Where's the algorithm? Um, it essentially produces the oiler walk backwards. It produces the last edge first. And here I I reversed it by putting the first edge I find at the end of an array. If you don't want to do that, if you want to output the walk in the correct order, you can just work on the reverse graph where you're using all the edges backwards and then the reversal is not needed. I find this confusing but uh it also works of course.

30:04 Initially all the edges are black and I'm going to encode what the algorithm does with uh two more colors. The red edge of each vertex is the first time it's reached in the traversal and then it might be reached many more times with other edges and these I'm going to also mark as visited u but in green and then the order among the green edges is irrelevant and the second important thing whenever I backtrack I have to prefer green edges over the red one let's see in the example and I think you'll you'll see how this maps to the original algorithm Uh, we're starting at this vertex. Now, remember, I'm doing everything backwards. So, I have to go to six. And this becomes a red edge because it's the first edge that reached six ever in its life.

30:54 Then I'm continue from six to five. This becomes another red edge for five. Now, here I have a choice, right? There's two incoming edges, but I'm going to two. And then from the two again, I could have gone back or I could go to the one. Let's go to the one. Vertex one I've seen before is the start vertex. So this edge is not red. Um the red edges always form a tree, the DFS tree in the standard sense. Um but now I have this other edge. I mark it as green as visited. Now at vertex one, it's the first time I have to backtrack. There's no more edges I can explore. Now this is where this green uh preference rule comes in. I have to follow the green edges. So backtracking usually takes the edges backwards, but here because we're working reverse, we're actually taking the green as the first. And so this will become the first edge in the final oiler walk. Okay?

31:49 And so we mark it as such. It gets the label one just uh for making the example easier to follow. And we mark it as dashed because this edge is essentially gone. We don't look at it ever again. Uh then we continue. We're back to vertex 2. Are there more edges? I could have followed incoming edges. Yes, I can go back to vertex 5. This is a green edge because I've been at five before. Then here I can continue up and up and up. So here another incoming edge, but I've seen two before. So that's going to get a green again. Now I've seen all the edges. So the rest will be all backtracking. And all I have to make sure is that I prefer green edges if there are any. Here there is a green edge to backtrack over. Here there's no choice. So I just do whatever. Um follow the edge. Here I'm going up for the green first and then red red.

32:43 Notice that there will always just be one outgoing red edge by construction because it's the edge from the DFS parent. There's only one ever. But green edges, I may have many more. In the tiny example, I never really had two green edges, but um to make it fit on the slide. And then I have my oilian walk here in this graph. Okay. Now, uh why did this reduce the space? This is from the colors. At least you can say I use I'm using m bits or order m bits to store the colors. But you actually don't need that. uh you can store everything implicitly um under the assumption that the input graph is represented in a way where it can iterate over incoming and outgoing edges. You can debate this but many many representations of graphs do this and many libraries implementing graphs actually explicitly duplicate directed edges because they want to support this traversal in both directions anyways. So I'm I'm just assuming that this doesn't cost me extra space.

33:54 But yeah, I can I can justify that. Um and as I said, every vertex has at most one outgoing red edge by construction. It's a DFS tree. So this one I'm going to store explicitly. I'm going to keep for every vertex an entry that says which was the red edge that I reached it via. And then the other edges that are visited uh which means they're red or green. I don't even have to store this at all because I'm like in the previous code. I'm just storing a pointer an iterator in the adjacy list of a vertex. And so all the all the edges that are after the iterator are not visited yet. All the edges before have been visited. And I never really have to ask is an edge visited or not. I just follow the iterator. So this costs essentially nothing extra in space. So all I need to store is uh the red edges and the iterator for the vertex and whether I've visited a vertex before. So these all all things that just take uh linear space and the number of vertices.

35:09 Maybe this might be good to say at the end. Okay. Um, give me one minute to wrap this up. Uh, if we look at where we came from in genomes, there's a few things that we haven't discussed um in in detail. uh the orarian walk might not exist because uh for some reason uh I have I have errors in my input or if it's not unique we may disqualify the graph as having a single walk. um more likely is that uh my input is just not a perfect k mirror composition and then I still want to do something and uh luckily I can still construct construct this deu graph and at least find what binformaticians then call context so is sequences of contiguous ranges of the genome that you can um that you can sequence even though you don't get the full genome in one go you can get larger larger fragments which are much much larger than the initial reads. And it seems that this is something people can then manually work with much better.

36:24 Um and we can get this from essentially looking at in-n-out degrees. Um this is just a summary of some challenges in practice. Um these are essentially for your interests. Um, so we don't really have a single string. It's broken up into chromosomes. So we certainly want to do one chromosome at a time at most. Most of us have two. So us in in sense of animals or ukareots maybe more generally we usually get chromosomes from the two parents and they're very similar. So this is actually horrible for sequencing. So the human reference genome is actually just one of the two chromosome pairs. And this seems to be uh acceptable. We might have gaps in the coverage of reads. That's a bigger problem. Um and we may have errors in the in the sequencing. And then another beautiful thing is the double strandedness. So we get essentially not just one text but the text and its reverse complement. And we'll get kr from both strands mixed. So that's that's another little puzzle to um to handle that part.

37:43 Uh there's a few things you can do. Imperfect coverage is relatively easy to address. If there's some starting positions you don't have, but you have long enough reads, you break them into shorter curs and then artificially you get much more reads than you started with. So this I thought was was really hard to imagine that the molecule is really broken up at every single possible base. But that's not at all an issue. You could just simulate that by making the reads a bit shorter and using all the possible starting positions within the reads. And the chromosome pairs also seem harmless. You can in the lab make sure you only sequence one. Uh the more in the more problematic thing seems to be the read errors and the de bruin graph completely dismisses if the overlap is not perfect. It completely says no edge go away. Uh and that seems to be the major trouble which is why people in practice tend to use the overlap graph because there uh it's much it's much easier to incorporate inexact overlaps give them a weight and then compute essentially a traveling salesperson tour of minimal cost. Uh where whereas the de bruin graphs with the lovely oiler walk theory seem less amunable to this generalization.

39:05 Okay, good. Uh that's that's all I wanted to say about uh the genome sequencing. Uh I think it's a it's a cute example of how the first obvious formalization of the problem runs into a hard problem and then if you look at it differently again it's actually a attractable problem.

39:27 Maybe the lesson there should be just because something looks np hard in the first attempt doesn't necessarily rule out an efficient algorithm. If you frame the problem differently,