Needleman-Wunsch global sequence alignment

[EDIT 01-01-2010] this is now available at bitbucket.
I've written a simple, fast, python version of Needleman-Wunsch as I couldn't find one to use. It uses Cython and specifically cython-numpy goodness. It's easy-installable as:
 sudo easy_install -UZ http://bpbio.googlecode.com/svn/trunk/nwalign/

or via svn from:
svn co http://bpbio.googlecode.com/svn/trunk/nwalign/

it will put an executable 'nwalign' into /usr/bin/ which when run will give this message:

Usage: 
    nwalign [options] seq1 seq2 
    

Options:
  -h, --help           show this help message and exit
  --gap=GAP            gap extend penalty (must be integer <= 0)
  --gap_init=GAP_INIT  gap start penalty (must be integer <= 0)
  --match=MATCH        match score (must be integer > 0)
  --mismatch=MISMATCH  gap penalty (must be integer < 0)
  --matrix=MATRIX      scoring matrix in ncbi/data/ format,
                       if not specificied, match/mismatch are used

where the matrix is optional but can be the full path specifying the transition score from row to column as here.
If the matrix is specified, match and mismatch are not used. If the matrix is not specified, match and mismatch are used. a typical run looks like:

$ nwalign alphabet alpet
alphabet
alp---et

$ nwalign --matrix /usr/share/ncbi/data/BLOSUM62 EEAEE EEEEG
EEAEE-
EE-EEG
And usage from a python script can be seen in test.py.

I wrote this for a colleague who was using a perl script. This is over 100 times faster than the perl script for long sequences--and the perl script had the BLOSUM62 matrix hard-coded as a hash. There are a couple places that could still be sped up, but I think the improvement will be relatively small. This kind of script is perfect for the numpy-cython mix as I can just access the contents of the 2d array at c-speed without having to do pointer arithmetic. If there's any obvious optimizations I missed, I'd be glad to know them.

Comments

Unknown said…
What is wrong with EMBOSS needle?
5:49 AM
brentp said…
@diffusing thoughts: the only thing wrong with it is that i didnt know it existed, and searches in google for various combinations of needleman-wunsch, global sequence alignment ..., etc. didn't find it. thanks for the pointer.
here's the link for anyone else interested.
6:57 AM
Michael said…
EMBOSS works well, but you have to a) write the sequences to files and b) do a process call/fork and c) read back the alignment/scores. I haven't tried this app yet, but it's a general problem (and Biopython only provides wrappers to EMBOSS so the problem remains)
8:58 PM
brentp said…
@Michael, that's good to know. let me know if you have any recommendations to make it more useful. i've since uploaded it to pypi
10:42 AM
Post a Comment

Popular posts from this blog

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

levenshtein in cython

EDIT(2): fix markup for <char *> casts ... fix malloc (see comments. thanks Bao). NOTE: using a kwarg for limit slows things down. setting that to a required arg and using calloc for m2 speed things up to nearly as fast as the pylevenshtein. Well, it seems to be popular to code up the levenshtein . I actually have a use for this and wanted to practice some Cython , so I've written a version. I used bearophile's recipe , wikibooks , this (from k4st on reddit) and this for reference. It follows bearophile's code closely, using only O(m) space instead of O(mn). cdef extern from "stdlib.h": ctypedef unsigned int size_t size_t strlen(char *s) void *malloc(size_t size) void free(void *ptr) int strcmp(char *a, char *b) cdef inline size_t imin(int a, int b, int c): if a if c return c return a if c return c return b cpdef int levenshtein(char *a, char *b, int limit=100): cdef int m = strlen(a),...
Read more

synteny mapping

Image
Living in Synteny I've been working on automating synteny mapping between any pairs of genomes. Synteny is where there's a stretch of DNA or genes in some order on chromosomeA of organismX and due to a shared evolutionary history, you can find a similar stretch of genes in order on chromosomeB of organismY. Often there are small losses and inversions, but between closely related organisms like man and mouse , there's still a lot of synteny. Plants can undergo polyploidy , following which, a species can have 2 entire copies of its genome. Over time, much of the duplicated cruft is lost, and the homologous chromosomes diverge, but if the divergence is not too great, it's still possible (actually common) to find synteny within the genome of a single organism--as well as between organisms. I've written my own algorithm to find synteny which uses python sets, and numpy array slicing to do the heavy lifting. It is quite clever [wink]. And it _almost_ works but it'...
Read more