We would love open source some of the work we did - there are a few edge cases to still work out with deprecated_column and renamed_to before I’d be comfortable doing that, but definitely agree that may be generally useful.
You're right - conceptually the CRISPR search problem and DNA sequence alignment are related. In both, you're looking for place where two (or more) sequences are very similar. I would say there are two major differences.
The first is in the goal of the search. Typically, alignment tools try to find the best positional alignment for two (or more) sequences. The CRISPR search problem is to find every possible match above some similarity threshold.
There are also a few constraints on the CRISPR search problem that allow us to make this much faster than a general DNA sequence alignment tool:
1) We know that that guide sizes tend to be very small (~20bp)
2) Part of the guide must match exactly (the PAM site), allowing us to restrict our search even further.
3) We don't need to worry about insertions or deletions in our search.
Using those three constraints, we can do this search a lot faster than a more general DNA alignment tool!
We actually do this in our "expensive" check. The reason that we use the 2 arrays is because in 2 array lookups, we can check for matches for all 200 guides. I probably could've made that clearer in the post - the array contains matches for ALL of the guides.
A lot of these are actually pretty easy to spot with tools like gcov. To determine false positives, we can just look at how many times each "if" statement was hit, and compare those to the final count.
Thanks for bringing up all of these alternatives. We definitely would have preferred to use an existing solution to building our own.
Unfortunately, a lot of the existing software is not intended for the search we're trying to do, or does not perform well under these conditions. We did in fact experiment with some of them before building our own. Bowtie, for example, doesn't allow more than 3 mismatches, and is also intended for alignments where there are very few matches (close to 1).
Since we need to be able to support multiple genomes (see Josh's comment), the amount of RAM we need to run a particular set of alignments is relatively important. Things like BLAT (which seems also intended for > 25 bases) need to keep the entire genome index in memory. This means that we would need to spin up a lot more servers to handle parallel requests, especially with different genomes.
FWIW our search is only a couple hundred lines of C++, and does the search with very little memory requirements.
We messed around with bowtie - it seems like most of these are optimized for the number of alignments being small (i.e. close to 1). Unfortunately, the number of matches for a 20 base guide on the human genome is closer to 1000.
I left PAM sites out of the blog post (it's actually mentioned briefly in the footnotes), as it made the problem slightly more complicated.
The final algorithm actually keeps track of the last 20 bases + PAM length, and checks both the edit distance and PAM before deciding if something is a match. The Benchling CRISPR tool will do this for you :)
We use the scoring function published by Hsu et al[1], which most scientists seem to be using. This function takes into account both the number of mismatches and where they occur in the guide. There's a more readable version here: http://crispr.mit.edu/about .
We would love open source some of the work we did - there are a few edge cases to still work out with deprecated_column and renamed_to before I’d be comfortable doing that, but definitely agree that may be generally useful.