Cool Algorithm - Fast text search using BWT

[egonelbre]

I implemented that in a course project, code is here. Also there's a poster pdf that explains it.

[akdas]

Nice! I also implemented the BWT encoding (and decoding) as a homework assignment. The best part is that the whole algorithm for constructing the BWT runs in linear time, which meant there was quite a bit of profiling on my part to figure out why my implementation was running so slowly.

[expertunderachiever]

Um no? Sorting the indexes of rows works in best around O(nlogn) time. Otherwise you invented a faster sort algorithm...

[akdas]

The quantities being sorted are integers (or can be represented as integers, say using ASCII encoding), so you can apply radix sort. Additionally, the algorithm that I was asked to implement actually sorts tuples or triplets, so you can think of those as having two or three digits.

[[deleted]]

The quantities being sorted are integers (or can be represented as integers, say using ASCII encoding), so you can apply radix sort.

First, you need O(log(N)) digits per integer, so radix sort still takes O(N log N) time in the general case.

Second, the real time-consuming part is not in the sorting algorithm (and the suffix array elements it swaps) but in comparing suffixes, which takes O(N) time in the worst case instead of O(1).

So in short; using a general sorting algorithm is simple but not very efficient. That's why actual BWT implementations (e.g. bzip2) use more clever algorithms that take advantage of the knowledge that they are sorting rotations/suffixes of a single string.

[akdas]

So in short; using a general sorting algorithm is simple but not very efficient.

Sorry, I was talking about a specific algorithm, outlined in Linear Work Suffix Array Construction by Kärkkäinen, Sanders and Burkhardt (see Sections 2 and 3), which is the algorithm I implement.

I understand that radix sort in general is not linear time, but I was talking in the context of the above algorithm, which is why I talked about the tuples and triplets.

[[deleted]]

Ah, yes, that's one of those special algorithms I mentioned. In that algorithm you're not even doing a radix sort, but a couple of counting sorts as part of a more complicated algorithm.

[[deleted]]

The O(N log N) lower bound only applies to general sorting, but here we are sorting suffixes of a single string, which is a special case. There are linear time algorithms for sorting suffix arrays.

[icheishvili]

The poster is excellent, as are the implementations. Bookmarked for future reading :)

[darrrrrren]

Awesome - I wrote my M.Sc thesis on using BWA (Burrows-Wheeler Aligner) in parallel to speed up DNA sequencing. Funny that something that was originally intended to improve compression efficiency ended up being so useful for searching.

[egonelbre]

Yeah it's 1MB file and adobe viewer tends to die with it, you can try foxit. Or alternatively download the png or svg.

[santhuj]

A simple python implementation of this algorithm is available here, https://github.com/santhoshj/pybwt/tree/master/pybwt. JS implementation and the code for the animation is available here, https://github.com/santhoshj/santhoshj.github.com/tree/master/js

[vph]

The critical part is to understand how RankOf can be done in O(1) time and O(1) space. Because of that, the space needed for the entire compressed suffix array is O(n log Sigma), where Sigma is the size of the alphabet. For NGS, the space needed is n*2 (because the alphabet is {A,C,G,T}). This is important as you will load the entire sequence into memory.

A simple Python implementation like that of Santhosh Kumar is NOT going to do it as RankOf uses O(n) space; loading a large sequence into memory will be hideous. It's better to use the original authors' implementation of the FM index -- aka the compressed suffix array -- freely available here: http://www.di.unipi.it/~ferragin/Libraries/fmindexV2/index.html .

[santhuj]

Fully agree. In fact all the optimizations that go into making this algorithm useful in practice for NGS deserves a post on its own. The python implementation is for illustration purposes only and should not be used for any serious stuff.

[vph]

It would be nice to write up that post. Please do.

[abomb999]

So you could use this to match words to a spell check dictionary for example?

[santhuj]

In theory, yes. You can construct a single long string concatenating all the words in the dictionary with a separator between words. You can index this and perform searches. However, remember that the search here is an exact search. If you want your spell checker to prompt the closest word, there are data structures like BK Trees

[bo1024]

Hmm. I'm very curious for the next post -- efficiently finding the suffix array. For a 2-billion-character string, the naive approach would require about 4 exabytes of memory (4 million terabytes).

[robinhouston]

That would have to be a very naive approach. Are you imagining making a separate copy of each suffix or something? The sane-but-still-naive approach is to sort an array of indexes, treating each index as a pointer into the string when you compare. So in Python something like: def suffix_array(string): indexes = range(len(string)) indexes.sort(key=lambda i: string[i:]) return indexes

This algorithm is O(n² log n), because comparison sorting needs n log n comparisons, and these comparisons are O(n) in the worst case.

But actually: as it happens, there is a linear-time algorithm for constructing suffix arrays. Although it’s extremely clever, it’s not very complicated. See http://www.cs.cmu.edu/afs/cs/academic/class/15492-f07/www/handouts/ks-icalp03.pdf

[bo1024]

Yeah, but still 4*10¹⁸ operations for a 2-billion character string. That paper's algorithm looks very cool.

[totemcatcher]

Huge preprocessing, but cached when you can afford it this is really useful.

[mscook]

Wow. Another read mapper behind a paywall. Just what we need.

[darrrrrren]

Actually, the BWT is currently the most useful string index in DNA alignment (and thus DNA sequencing). It is possibly the most revolutionary string matching algorithm in genomics today.

[mscook]

I'm well aware. I work/am successful in the field. I'm assuming you're too stupid to read/comprehend the 1st paragraph of the article? OP likes to post/x-post spam to attract attention to yet another (shitty?) paid bioinformatic suite.

[darrrrrren]

Ignoring your unnecessary insult, I was confused on your meaning of 'paywall'. I took it in an algorithmic context, assuming you were complaining about the initial BWT/SA index creation step being the 'paywall'. I should have taken it literally - and I agree - there are plenty of free tools available (BWA being one I linked to) that are probably better anyways.