Genedex: query genomic features and sequence

Normally, I don't write libraries, I figure smarter people than I should do such things, and I should just use them. But, I got tired enough of writing one-off scripts for genomic feature manipulation-- find the upstream, downstream neighbors and get the sequence -- and I saw enough of the pieces coming together that I decided to build it. I'd complained before about how rtree didn't support 1D indicies. Not only is this not a problem, it's beneficial. Genomic features should have strand information, so that's the 2nd dimension. Then rtree does containment queries, so it's simple to find only the features on a given strand. I realized this about the same time that the docstring for numpy's memmap went from 0 lines to about 100, and it was enhanced to take a filehandle, not just a filename. This means you can send in a start position and a shape to the numpy.memmap constuctor and it can create a numpy array of only that chunk. This means that it's possible to slice an unaltered fasta file using the numpy array syntax. That's very good.

So, if you put those 2 simple things together, you have the start of something powerful. That's what I did. Then I gave it a crappy name: Genedex (Gendex was taken) and slapped it into googlecode. Check it out: http://code.google.com/p/genedex/. My only design goal was to keep it as simple as possible. If the amount of features is under-whelming, that's good.

TDD


Also, I generally do TDD very half-ass, with asserts and maybe a couple doctests. However, I recently made fairly substantial changes to the SQLite datasource in featureserver, and wrote this set of doctests while doing so. It works! and I've been using it. So, I did what featuresever (presumably crschmidt) devs did and copied the setup for the shapely doctests. It's pretty useful for design, i'd just write out the code for how I wanted the API to look and then implement. The only thing is, for doctests, the way they're used (at least by me) is to copy the output from executing the code into the doctest. So, if your code is wrong to start with, you just copy the wrong answer into the doctest and it's broken but the tests pass. But, at least it's good for regressions, and I just had to remember not to blindly trust the output. That's true for all testing, but especially so for doctests.

So, there's now more tests than code. But, since it's mostly just tie-ing together pieces that do the real work, it's not much code. Doc-tests are also nice because (as the name suggests) they double as documentation. So, here's the genedex documentation:
http://genedex.googlecode.com/svn/trunk/doc/readme.html
It's pretty! It gets colored by pygments, using this script. The only major thing I'd like to add to the library is a plotting class using matplotlib. Then other smaller tasks like a method that takes 2 features and returns the sequence between them.
Any fixes, enhancements, ridicule, etc. will be greeted with commit access.

Comments

Post a Comment

Popular posts from this blog

filtering paired end reads (high throughput sequencing)

NOTE: I don't recommend using this code. It is not supported and currently does not work for some sets of reads. If you use it, be prepared to fix it. I wrote last time about a pipeline for high-throughput sequence data. In it, I mentioned that the fastx toolkit works well for filtering but does not handle paired end reads. The problem is that you can filter each end (file) of reads independently, but most aligners expect that the nth record in file 1 will be the pair of the nth record in file 2. That may not be the case if one end of the pair is completely removed while the other remains. At the end of this post is the code for a simple python script that clips an adaptor sequences and trims low-quality bases from paired end reads. It simply calls fastx toolkit (which is assumed to be on your path). It uses fastx_clipper if an adaptor sequence is specified and then pipes the output to fastq_quality_trimmer for each file then loops through the filtered output and keeps onl
Read more

python interval tree

EDIT: added a couple points inline. I'm obsessed with trees lately -- of the CS variety, not the plant variety. Although we are studying poplar, so I'll be using trees to study trees . I'd tried a couple times to implement an interval tree from scratch following the Wikipedia entry . Today I more or less did that in python. It's the simplest possible form. There's no insertion (though that's trivial to add), it just takes a list of 'things' with start and stop attributes and creates a tree with a .find() method. The wikipedia entry that baffled me was about storing 2 copies of each node's intervals --one sorted by start and the other by stop. I didn't do that as I think in most cases it won't improve lookup time. I figure if you have 1 million elements and a tree of depth 16, then you have on average 15 intervals per node (actually fewer since there are the non-leaf nodes). So I just brute force each of those nodes and move to the next. I th
Read more

needleman-wunsch global sequence alignment -- updates and optimizations to nwalign

I've written previously about nwalign , a python package I wrote to do fast (thanks to cython) global sequence alignment. Over break, I spent some time fixing some bugs and improving performance. Ack It's actually nice to get bug reports for a relatively obscure bit of software like this as it shows it's getting used. Thanks especially to R. Christen for his patience in (repeatedly) showing me places where nwalign was not doing the right thing. Bugs Placement of Gaps Some of the "Bugs" were actually just "unspecified behavior". For example, given input text of "AGEBAMAM" and "AGEBAM", the alignments: AGEBAMAM AGEBAM-- and AGEBAMAM AGEB--AM Have the same score. However, the first version is generally more, um, polite. Previously, the behavior was deterministic, but depended on the length and order of the input sequence. As of version 0.3, nwalign will always put matches together when given a choice between 2 (or more) paths through t
Read more