BHLog

Roche’s 454 sequencing technology can produce biological sequence data on a scale that exceeds traditional Sanger sequencing by orders of magnitude.   Due to the fundamentally different method used to generate the sequences, we would like to investigate the raw data and see if we can quantify — and maybe also reduce the number or severity of — errors.  This means reading the binary SFF format. Below, we’ll dissect the SFF format, and describe a Haskell implementation.

Read more…

I’ve been spending the last couple of weeks working on an indexing scheme for sequences, using Bryan O’Sullivan’s Bloom filters.  Now, it turned out that when Bryan tested out the code, he found a curious problem:  Apparently, the indexing stage scaled quadratically with sequence length.  This wouldn’t have been so strange, were it not for the fact that I saw the expected linear time usage when I ran the code.  Some more digging about revealed that my laptop also showed quadratic scaling.  Profiling showed the  culprit to be a simple pipeline-style function:

swords s = take (fromIntegral (seqlength s)+1-k) . map (B.take k) . B.tails . B.map toUpper $ seqdata s

Here, seqdata returns a lazy bytestring, which is also what’s hiding behind the B. qualifier.  Basically, this just builds the list of all lenght-k substrings.  This should, in my opinion, stream nicely, and run in constant space and linear time.  What on earth makes this quadratic?  On only some systems?   Time to dissect the systems in question:

Read more…

It was just brought to my attention that people have started to use a new file format for sequences.  This format, called ‘FastQ’ combines both the sequence data itself and the quality data in one file.  That’s a nice idea, and I implemented support for it, tests, docs and all, in the bio library.  Runs fast, too.  Basically, the format is a sequence of records, each one similar to this:
Read more…

In my copious spare time, I’ve started putting together a tutorial for using the biohaskell library.  It’ll probably take some time — anything from a while to an eternity — until it’s complete, but I thought I’d follow the adage of “release early, release often” in the hope that the intermediate product may prove useful to somebody, somewhere.  It shouldn’t really require any prior knowledge of Haskell nor of biology.  Anyway, it’s here, please take a look, and tell me what you think!

In other news, inspired by an article on programmer productivity at LWN, I ran darcs-graph on my code to see how I do.  I guess I consider myself an about average programmer, and it looks like I can average about five commits a day when I’m working on a project.  I’m occasionally touching 20 commits, but that’s probably a built-up backlog of patches.  Let’s see where this puts me:

Read more…

Bloom filters is apparently a relatively old technology, dating from the 1970s or so, but it has somehow escaped my radar until Bryan O’Sullivan posted a message to the haskell mailing list announcing a high-performance implementation in Haskell, perhaps to support a chapter in the upcoming book.  You can read all about Bloom filters on Wikipedia, but the executive summary of it is that it is a structure similar to Data.Set.  Except that it is probabilistic, and may occasionally claim a value is a member when it’s not.  On the positive side, the Bloom filter is very fast, and speed is independent on the size — in other words, lookup and insert is O(1) where Data.Set is O(log n).

Comparing sequences to find similarity is a common occurrence in bioinformatics.  For instance, one might want to know where a certain gene is located in the chromosome, or which sequence fragments are similar enough to originate from the same gene. To speed up searches, it is common to index sequences in questions as overlapping, substrings (k-tuples, q-grams).  This index seems like an obvious target for Bloom filters — large data, time critical, some false positives anyway — but for some reason, there is almost no such applications that use them. Until now.

Read more…

Just some quick notes:

Hackage submissions updated

There seems to have been problems with some of the bioinformatics applications on Hackage, thanks to Don S. for pointing it out. That should be fixed now by new uploads, but I’m still waiting for the automatic builds to register results. An, since you ask, it was all my fault for being sloppy with version dependencies. It’d all work with a recent biolib. Speaking of which,

A home page for the bioinformatics library

I’ve finally updated the static home page for the library, it can be admired (especially if you remember what Oscar Wilde had to say about that) here.

I’ve discovered HPC

No, not high-performance computing (which sucks anyway), but GHC’s new ability to do coverage profiling.  Adding it to the default testing procedure was just a couple of extra options, the results can be admired by browsing the HTML files in the darcs repository.  This is very cool, almost too easy, and looks pretty good. The only downside is that it exposes my sloppy attitude to testing.

My current development project is an EST pipeline.  For various reasons, it is implemented in shell — bash, to be exact.  In other words, the pipeline is a script, or rather a set of scripts, that will tie together the various stages: masking, clustering, assembly, and annotation.

As in any program, there are many occasions where you want to effect some particular change during some part of the program.  The archetypical example is allocation of local variables. After allocation, the variables are then available to the program until they run out of scope, they then get deallocated automatically.  The technique can be generalized beyond this.  For instance, you (or rather I) may want to set a $STAGE variable that indicates the current processing stage, and which should be unset when the stage has finished executing.  Or you may want to run some processing in a different directory, in which case you really want to remember to return to the previous directory when you finish.  The purpose of bracketing is to wrap a section of code with an initial part to be run in advance, and a final part to be run afterwards.

Read more…

As a consequence of IWC policies, the Institute of Marine Research is required to store genetic identification of each minke whale that is hunted. This of course means that people will come to me for help in bridging the gap between test tubes and the databases by providing some analysis tools that can extract the information from the data.

Well, I’m not complaining, work security and all that, but it does leave me with a slight problem. I’m a Linux user, and my experience with non-Unix OSes after I passed on my Amiga 500 to my little brother is rather limited. (It’s not that I won’t fix your PC, rather, it’s quite likely that I can’t.) But the users I support in this particular case are largely confined by the walls of Redmond, and not terribly enthusiastic about broadening their horizons in that regard. Here’s one way to solve it. Read more…

It seems to be optimization week on the haskell café mailing list. Inspired by a combination of Don Stewart’s blog post about how to optimize for speed and the sorry performance of my xml2x program, I thought this would be good time to see if things could be improved. In its current state, xml2x takes a few hours to process a few days worth of BLASTX output, so it’s far from critical, but faster is always better, and reading Don’s post, the intermediate output from the compiler, a.k.a. ghc core didn’t really look so scary. Read more…

…to write a simple, but useful and efficient bioinformatics program? Here’s how to build a simple tool to extract a clustering from an ACE file, using functionality in the bioinformatics library. It is about three lines of actual code, it is fast, and it is efficient.

Read more…