A new online viewer for multiple alignments

For my work on the alignment of coiled-coil proteins, I needed an alignment viewer that could highlight coiled-coil domains, since they contain less phylogenetic signal than other parts of the protein. Adding this to JalView seemed very complicated, and not possible at all in other web-based viewers like MView. I therefore created amview (for "annotated multiple alignment viewer"), which looks like this in practice:

Shown is the MSA for spd-5 (and here is the entry in my coiled-coil orthologs database). Amino acids are colored according the ClustalW rules, coiled-coil residues are in a lighter color (and the "a" register of the heptad repeat is underlined). Below, two interaction domains are shown. At the very bottom, a small chart shows the degree of conservation across the whole alignment, which can be used to quickly scroll to the conserved regions. The cog on top hides two options: you can hide columns with too many gaps, and proteins that seem to be fragments.

I've only tested this thoroughly in Google Chrome, as I found other browsers to be too slow. Still, it's better than loading a Java applet, and even runs on iPhones/iPads etc.! The implementation relies on Django on the server side, and JavaScript / JQuery in the browser.

New repo: local NCBI taxonomy database

I added some new functionality to the taxonomy repository at Google Code, creating a fork at BitBucket. The existing Python package already makes it possible to create a local database containing the content of the NCBI taxonomy, which can then be queried for names, ranks, and lineages. I added functionality to create a Newick tree from a list of NCBI taxonomy identifiers.

A Newick parser for Python, supporting internal node labels

I just pushed a fork of Thomas Mailund's nice Newick parser for Python to bitbucket. I added support for labeled internal nodes, but probably partially broke support for bootstrap values.

>>> from newick import parse_tree
>>> t = parse_tree("((Human,Chimp)Primate,(Mouse,Rat)Rodent)Supraprimates;")
>>> print t
(('Human', 'Chimp')Primate, ('Mouse', 'Rat')Rodent)Supraprimates
>>> print t.identifier
Supraprimates

Using Makefiles for jobs that run on a cluster

Makefiles are great. While you work on a project, they make it convenient to run the necessary scripts. When you come back to the project half a year later, you don't have to dig in your brain how the scripts fit together—it's all there. (More on make, and related advice.)

However, often in bioinformatics computational tasks are too big for a single CPU, so jobs are submitted to a cluster. Then, the Makefile doesn't help you much: It can't detect that jobs are running on the cluster. There is qmake, but it only works if all your parallel jobs are specified in the Makefile. I usually write my parallel scripts in a way that they can submit as many instances as necessary of themselves to the cluster via qsub.

Therefore, I went ahead and wrote a small Python wrapper script that runs the job submission script and sniffs the job ids from the output of qsub. It then waits and monitors these jobs until they are all done. Then, the execution of the Makefile can continue.

Here's an example of how to invoke the wrapper script from the Makefile:

pubchem_compound_inchi.tsv.gz:
   ~/src/misc/qwrap.py ${SRC_DIR}/inchikeys.py
   cat ../inchikey/* | gzip > pubchem_compound_inchi.tsv.gz

You can download the code (released under a BSD License, adapted to SGE). I hope it's useful!

Addendum: Hunting around in the SGE documentation I found the "-sync" option, which, together with job arrays, probably provides the same functionality but also checks the exit status of the jobs.

Sending Growl notifications from Python scripts

Working in bioinformatics can be seen as an infinite loop of: think, write a script, run script, analyze data. While I work on a Mac, most of the scripts run on a Linux server, and it would be nice to know when a script is done so that I can look at the data. In order to be notified when a script finishes, I now use Growl (see picture).

I downloaded netgrowl.py and wrote a quick and dirty wrapper package around it:

#!/usr/bin/env python

from netgrowl import *
import sys

def growlNotify(title = "Script Finished", message = ""):
    
    addr = ("10.1.104.26", GROWL_UDP_PORT)
    s = socket(AF_INET,SOCK_DGRAM)
    # 
    # p = GrowlRegistrationPacket(application="Network Demo", password="?")
    # p.addNotification("Script Finished", enabled=True)
    # 
    # s.sendto(p.payload(), addr)

    if not message:
        message = sys.argv[0]

    p = GrowlNotificationPacket(application="Network Demo",
        notification="Script Finished", title=title,
        description=message, priority=1,
        sticky=True, password="?")
    s.sendto(p.payload(),addr)
    s.close()

if __name__ == '__main__':
    growlNotify()

The registration is hidden in a comment, you only need to do that the very first time. So, in my scripts, I just insert the following right before the end of the script (using a Textmate snippet to save typing).

import growlnotify
growlnotify.growlNotify()

Use case for decorators

Intro

Decorators were introduced in Python 2.4 but remain a somewhat obscure feature of the language. A decorator essentially is a function that you apply to another function, replacing it with the return value of the decorator function:

def d(f):
...

@d
def g():
...

This is equivalent to:

def d(f):
...

def g():
...

g = d(g)

More background here.

The use case

As part of my research project, I’m testing several ways to explain side effects of drugs. So I made a class that loads all the relevant data into memory and that contains methods to fit the data. (I thought about creating subclasses replacing the “fit” method of the base class. However, the data that’s being loaded into memory is the same for all instances of the subclasses, so I don’t want to reload it every time. Thus, I put multiple methods into the base class.)
At first, I used introspection to identify the relevant functions by name (using dir) and then calling them. This worked fine, but it doesn’t make it particularly easy to disable some of the functions. Enter decorators: I wrote a decorator that takes a function, puts it into a list (which is a global variable) and returns the function unchanged. To call every fitting function, I can just iterate over the list. If I want to disable a function, I just comment out one line. Yay!



Another neat trick is that the docstring of each function is accessible as an attribute, so that each fit comes together with a description.