Video 6-3: Generalized suffix trees

0:02 Now let's try to come back to the original problem that I motivated this section with the uh longest common substring between a collection of strings and we'll see a key uh tool that we can apply to suffix trees in general here that is very very useful in in all sorts of problems.

0:25 what we uh seem to need. If you compare what we currently looked at, uh the longest repeated substring of a single string and the longest common substring of a collection of texts, it feels like we're solving almost the same task, just that here it's several strings and here is one string. So, we somehow uh need a way to compare across strings. So you could build the suffix tree for each of the texts, but that doesn't seem to help you find matches between the texts. Somehow we need a way to get that all in in a single uh data structure. Now the suffix tree is strictly defined for a single string. But what we can do is concatenate the strings. And um you want to be a little careful so that you don't have mixtures um a suffix of one string and a prefix of another string. um these we are not interested in. So we put these extra end of string markers and they're index they're index so that they're all all different um between the the different texts. The order doesn't matter it's arbitrary. So we build one long text by concatenating all the texts and separating them with a unique symbol. And then for that we built the suffix tree. So that means whenever you start in in some some text. Then there's a dollar, there's another text, another dollar, and so on. If you say we're looking at this suffix here, then really once you see the first dollar symbol, you can you can stop because you know from then on uh this will never match anywhere else. So this will always be one long edge towards a leaf. Right? As soon as you see this, the path needs to be unique because this dollar symbol occurs nowhere else. It cannot occur again and uh there can can be no sharing beyond this point.

2:28 So uh we will have some below below this initial part that is part of of a single TI. uh this will branch off into some some leaf and that will always be the entire rest of the string.

2:47 What that means is for x for every suffix and every for every suffix of every text we'll have exactly one leaf and we can pretend that they're labeled like this. uh they should really be labeled with I as the starting index in the entire uh they're originally labeled with starting indices in this concatenated string but you can just recmp compute what uh the text is where it starts and the position in that text. So we'll ignore how to compute these you just subtract the cumulative sums.

3:23 Here's an example for these three uh texts um all over the the alphabet ABC. Um and what you notice is there's a whole bunch of leaves, but uh I've always stopped after the first dollar symbol. So this would normally where if you directly apply the construction, uh this first one would continue with the whole second text and dollar two and the whole third text and dollar three. But we don't care about this. We all all we want is the suffix of the first text.

4:00 Other than that, it's using the same construction. You're just building it for this concatenated text and superficially I'm rewrite the leaves and remember the leaves are not stored as strings anyways. These are just starting indices. So, uh we're just reinterpreting these numbers.

4:21 So, what can we do with this? Uh, one thing we can do is the longest common substring problem. That was this original one where Don can at some point conjectured it might not have a linear time solution. Uh, but it does. And here it is. You compute this generalized suffix tree that we just introduced. And then for every internal node, you write down what kind of leaves are below this node. Namely, kind means what substring, what subtext, what text do they correspond to? And there might be some internal nodes that only have one type of leaf. They all belong to the same text. Or it could be any subset of the strings. So we do that bottom up. For all the leaves, we just write J because that's the the single string that this leaf corresponds to. And then bottom up when we've done all the children, we take the parent and we just compute the union of the two children.

5:14 Uh we also need the string defs as before. And then the longest common substring that is common to all k strings is this the internal vertex of of largest string depth that is labeled with all strings and that's it. Uh before we discuss the running time let's see this um example.

5:37 So, uh, here's, um, the same suffix generalized suffix tree from the slide before. And here's this annotation with taking the leaf types. So, uh, this leaf dollar two here is of type two. This is of type three. This is of type one. So, this node already has all strings below it. That's why it has one, two, and three. whereas some others might have just a subset. So among all those that um have all three that will be some top part of the tree we find one that has the largest string depth. In this case this is uniquely this one. It has a three. Everything else is smaller. So BCA is a common substring common to all three. And there it is. Right.

6:34 And it's also the longest such shared string. Otherwise, we would have another node. Uh in the code, um most of these steps just take linear time in the total length of all strings. In particular, computing the generalized suffix tree is computing a suffix tree on the string of all concatenations. So that takes takes the sum of all the lengths. Um and these are also uh simple traversals that just look at every node in the suffix tree once and the number of nodes in the suffix tree is at most n. Now this uh second step um if you do it uh just as I described it then at each node we might have to deal with sets of size up to k. So you probably want to represent them maybe as a sorted sorted list or maybe some bit vector. Either way, um you might need to spend uh as much time as there are different strings there per node. So in particular, if the number of strings is small and they're just very long, then this is overall linear.

7:44 If K is bigger, then in general, the time is K* N in this implementation. Okay, so that's essentially this original problem from Duncan solved and the main ingredient that we needed is the suffix tree.

8:14 This construction that we did uh doing the generalized suffix tree and then annotating nodes by what kind of leaves are below them is very versatile. You can solve all kinds of similar problems uh using the same approach. And one I put this slide in green because it's it's really more u motivated from lab techniques is to check whether your sample of reads from some some DNA sample is contaminated with a DNA sequence that you're you're not really interested in but which was accidentally dropped into the tube in the lab. So um before we formally model this, let me motivate a bit that this indeed can happen. Uh here's an article from uh this journal science.

9:04 Where are we? uh which is about some protein sequences from some uh some mammoth but more famously from a T-Rex. So 68 million years old. That's what it says here. Uh and it's it's the story of uh a scientist um was a student of one of the most famous paleontologists uh of of our time worked a lot with uh dinosaur bones bones and she was uh she was hoping that uh a bit of bone from this T-Rex skeleton that they found and that she accidentally dissolved a bit longer in some demineralization. solution than intended. There was some uh some substance left that looked a bit like blood vessels. And so she thought, "Okay, cool. There might actually be uh nonbone substances, proteins from this T-Rex left after all those many years, 68 million years. So uh they sent it to through mass spectrometry and applied some uh reconstruction techniques.

10:33 uh without getting into the details of how that worked um they tried to determine the sequence of the peptides of the proteins in this uh sample.

10:43 Now this caused a lot of spur. You may know nature and science as the the top journals in in the sciences natural sciences. So to get a paper there you have to show something surprising and this certainly made some some uh uh news at the time. Um but then there were several reanalyses of the same data and also uh suggestions of what might have happened. I think it's not completely clear. Um but here they at least claim that it might well be that the sample contained common laboratory contaminants like soil bacteria, birdlike hemoglobin and collagen. And in particular at the time, this first paper detected a surprising similarity between the proteins in T-Rex and chicken, which was taken as the first molecular justification that indeed birds are the successors of of dinosaurs, which is is otherwise not contested, but it was the first molecular such proof. But it turns out that at least the the evidence is very shaky.

11:56 And so it's it's kind of a real problem. Uh you can easily imagine how uh some uh well you're we're full of protein, right? Whatever exactly um happens in the lab, even if it's just you working on it, that there's some contaminants going in there is is very possible. So as the second article said, there's some likely contaminants that often end up in samples like this. I guess if you dig up dinosaur bones then they're often in soil. So having soil bacteria peptides in there is something you could expect. [snorts] So suppose we know some typical contamination sources and now we check our new reads. In this case it's formulated for DNA not protein but the principle for the abstraction is the same. So we have some some reads. These come from some uh DNA sequencing run and we have some known sources from a database of what typical I don't know in this case maybe chicken protein uh looks like [snorts] and we'd like to find for all those reads only those that don't have a long common substring with any of the contaminants because if they do then maybe we should discard those reads as being a contamination source but we don't want to. So the whole goal here is not to throw the entire sample away, but just keep the reads that are likely new and not from any of these known sources. And we're we're saying here common substring because it's it's enough for us if there's a a decent overlap. We know that the reads are not uh precisely read and we may have only a snippet of some longer gene. So it's enough if there's a long common substring that we declare this a a potentially a a contaminated read. So how how would we do this? Um it looks first of all a little different from the previous problem but you can cast it in the same framework.

14:02 Uh here we built the generalized suffix trees of all those sequences all the contaminants and all the reads. Uh in practice this might be a problem to do this in one shot but we we'll ignore this for now. Then we uh we look at one of those reads PR. This has a common substring with some specific so fixing one one read and one contaminant sequence. This has a sub a common substring of length L. Well, exactly if there's um think about it. If the this is our PR and here's our contaminant uh CJ, there must be if there's a shared subsequence, right? We don't know where it starts, but this means there's a suffix starting here and suffix starting here of the two sequences that share a common path in the suffix tree. So they're the same up to a point and then they split. So these two points here, this corresponds to an internal node in the suffix tree where they branch off. So we can find this by simply looking um for internal nodes that are at sufficient string depth. We only care about long shared uh substrings.

15:21 Um and they need to have any leaf coming from these two strings below. Then that's a a match. we have a long substring between the two. So if we mark all the um nodes in t where with uh it's just a binary decision is there any contamination leave below you then you get a a flag as a dangerous node. Uh, and then for any of these marked nodes, we check if they have a a high string depth above L, then any read that has a leaf below them needs to be removed from the sequence from the sample.

16:07 And so uh if you do this smartly, you can do it in linear time in the total length of all those sequences, which is uh is pretty surprising if you think about it from first principles because it it feels like you have to check every pair, but the the cleverness of suffix trees solve it for us.

16:31 As a last example of this kind, um we'll revisit a structure that we've seen um a while ago where we were talking about the the sequencing problem for genomes where we want to take all the reads and assemble them into one long string. Um there we were trying to do it with perfect overlaps so that you have k mirrors all of the same length and then there's exactly one character that we're missing here and missing here and these are all in perfect correspondence that was the setup at the time. If that is the case then you can use the de Bruyne graph solution we've shown there with oiler cycles and the whole problem is um is efficiently solvable but um with inexact reads we can't expect that this overlap is perfect and so uh one option that we have is instead of having this perfect overlap we have sequences of uh then we might as well have arbitrary lengths and we just find the longest overlap where I can put one sequence on top of the other and um if if they don't match too well then this is maybe just a short overlap by a few characters whereas if they're almost matching like this then probably they overlap by quite a bit that's the idea that means um in the end we'll build a graph of these sequences this Uh which has a weight that corresponds to the number of characters that we overlap here or uh number of of mismatched characters. Um depending whether you want to minimize or maximize So if you pick out this example here, say if it goes from here to there, then maybe this is the suffix we can find as a prefix here.

18:48 And then we're back to the TSP or or Hamiltonian path problem. We try to find a round trip that uses every uh read only once and where the total uh in this case the total number of mapped overlapped characters is maximized. Now that second step is in in theory a hard problem. uh what we'll focus on here is uh computing all those weights because we have hundreds of millions of reads in a usual uh full genome sequence uh run. If we do all pairs of reads and try to find these overlaps that's never going to finish. This is completely beyond it's out of the question. So how can we compute um all these overlaps? Formally we say we have a bunch of reads and we're interested in this matrix uh of overlaps where this is the length of the longest suffix of one read that uh overlaps with the prefix of the next read.

20:01 So in this formulation uh the output is already quadratic in size. So you would have to spend quatratic time. What you would want to do in practice is uh any any kind of entries that are below a certain length you just discard and don't make the edge there. Um for illustrating how we use suffix trees to do all this uh overlap computation. Um let's ignore this part. Uh we'll build the generalized suffix trees over all the reads as before. Now we have to find somehow um let me show the example first. I think this is a bit easier. Um here's the the same same suffix trees as before. And let's try to look at at this T2 this string as the um second string. Uh so I'm looking for prefixes of T2.

21:08 So there's these other sequences that uh come before and then there's T2. And um I'm interested in the longest uh prefix here that I can find as a suffix in the other strings. For example, if you look at BC AA, then this edge should have a a length two because there's this AA that we can match. And similarly for the others, how can we use the suffix tree to find this? Well, the suffix tree contains all the suffixes of the other strings. So, these are already in the tree. All I need to find is these prefixes. But the prefixes I get by traversing down this string and stopping at every position if you like. So we'll take this string maybe a a bca. I'll start at the root and go down a a and then bca. Then I'm at the leaf. And whenever at this journey there is one of those green marked leaves namely those that just have a dollar symbol that means there's some other suffix of some other string that is essentially done at this point apart from the end of string marker.

22:36 So uh uninterestingly this is true in the root I can always have the empty overlap. So, uh, there's all of these strings potentially present many times, and I have to find the deepest of them, the longest. So, I go down so far from for all these, uh, for all three, I've only found an overlap of zero. Amazing. I go down one character. Now, for two and three, I've now found a longer overlap, namely one character, the A.

23:06 If I go down one further step, I find oh for three it can actually do two characters and indeed I can also do two characters for two itself. Maybe is a bit uh amusing. No, it's not right. That that's not a terminal link. That's not green. I have to stick with the two here. I also get a self loop right in this in this way um of uh length one. We're not going to use these in the sequencing uh application, but you could say the algorithm just computes this for free. Uh and from the other string um clearly because it ends with a C, this is going to have an edge of weight zero or maybe the edge is not even there. So the algorithm has to do two things. It has to uh label the edges that have just a dollar symbol somehow specially and then it has to traverse every string from top to bottom and take note of what was the deepest occurrence of each other string on such a green edge.

24:19 That's what's spelled out here. These uh green edges are called terminal edges. They're the ones that start directly with a a dollar. symbol and um we have a a shared suffix prefix suffix of pi is is a prefix of pj. If that traversal of pj in t reaches such a terminal edge and we want to find the deepest of these. So to do it cleverly, you do a single traversal um with uh each of the of the reads and for all the other strings, you keep track of what was the deepest you've seen so far. And then when you uh when you reach a leaf, then you output all of these. And you can do you can even do a single traversal of the entire tree uh not just for any specific read, but you can do all of the traversals. um together. Then the overall time um is is linear plus the output size.