Python for Bioinformatics

Wednesday, November 17, 2010

simple text in matplotlib

This post is a "note to myself" about fonts and text using matplotlib. It can be hard to remember things like what is a FontProperties object and what kind of object is a font_dict, so here is a little example. It shows how to get (1) the default font with italics, (2) a default font in a particular "family"---sans-serif, (3) and (4) two named fonts and (5) a font loaded from a particular .ttf file. Notice that in (4) we do not observe the desired response to changing the 'weight'---I'm not sure why yet.

The second figure shows a change in the font using rc. It looks nice, takes several seconds to draw, and does not give the desired italics. They're all here as a reference. Perhaps they'll be of some use to you as well.

output:

:family=Arial:style=normal:variant=normal:weight=normal:stretch=normal:file=/Library/Fonts/Arial Italic.ttf:size=12.0

import matplotlib.pyplot as plt
import matplotlib.font_manager as fm

# can take a 'fontdict', Python dictionary
x,y = 0.1,0.1
fd = { 'fontsize':24, 'style':'italic',
       'va':'center', 'ha':'left' }
plt.text(x,y,'1. text ' + 'default',**fd)

fd['family'] = 'sans-serif'
y += 0.1
plt.text(x,y,'1. text ' + 'sans-serif',**fd)

L = ['Arial','Helvetica']
for i,n in enumerate(L):
    fd['fontname'] = n
    y += 0.1
    plt.text(x,y,str(i+2) + '. text ' + n,**fd)

fd['weight'] = 1000
y += 0.1
plt.text(x,y,'4. text--not properly bolded',**fd)

fp = fm.FontProperties(family='Arial',
    fname='/Library/Fonts/Arial Italic.ttf')
print fp
y += 0.1
plt.text(x,y,'5. text Arial--named',
        fontsize=24,
        fontproperties=fp)

ax = plt.axes()
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
plt.savefig('example.png')

from matplotlib import rc
import matplotlib.pyplot as plt

x,y = 0.1,0.1

rc('font',**{'family':'sans-serif',
             'sans-serif':['Helvetica'],
             'style':'italic',
             'size':24 })
rc('text', usetex=True)
y += 0.1
plt.text(x,y,'text using rc--not properly italic')


ax = plt.axes()
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
plt.savefig('example.png')

Simple server running Python script on OS X: part 4

Continuing with the project from last time and before (here, here, and here).

I got two more things working: a dropdown menu and sending a file.

It would be nice if the user could only send the file types we want to accept, but that's not the way it works, at least you can't disable the unwanted types in the file dialog the way you can from Cocoa. I tried to use an accept attribute, but then I found it isn't supported by any browser! (Lost track of where I saw this, it's not mentioned in the spec).

Two possibilities: javascript to check the file extension, or coding in the script. I wasn't able to find MIME type info, see the printout. Probably the best would be to use Python-magic or something like it.

I had some errors that seemed to be related to file permissions, but at the moment it is working with:

$ ls -al data.txt
-r--------@ 1 telliott_admin  staff  41 Nov 17 08:23 data.txt

I'm having trouble with output, so here is a screenshot of the first few lines:

form.html:

<Content-type: text/html>

<form name="input" enctype="multipart/form-data"
action="cgi-bin/script.py" method="post">

<select name="dropdown">
<option value="Math" selected>Math</option>
<option value="Physics">Physics</option>
</select>
<p></p>
<p>Please input a <code>.txt</code> or <code>.html</code> file:</p>
<input type="file" name="filename" />
<p></p>
</select>
<input type=submit value="Send" />
</form>
</html>

script.py

#!/usr/bin/python
import os, cgi
import subprocess
form = cgi.FieldStorage()
D = os.environ

s = '''Content-type: text/html

<head></head>'''

print s
print form.keys(),'<br>'
print form.getvalue('dropdown'),'<br>'
data = form.getvalue('filename'),'<br>'
print data[0],'<br>'
print len(data),'<br>'

print '<table border="1">'
for k in D:
    print '<tr>'
    print '<td>', k, '</td>',
    print '<td>', D[k],'</td>'
    print '</tr>'
print '</table></p></p>'

print 'CONTENT_TYPE'
print cgi.parse_header('CONTENT_TYPE')
cgi.print_environ_usage()
print '</body></html>'

Tuesday, November 16, 2010

Simple server running Python script on OS X: part 3

Continuing with the project from last time and before (here and here).

I got the form to work. You can use GET as long as there isn't too much data, otherwise, you should use POST. With GET it comes in the QUERY_STRING, as shown last time. The specs are here and some useful info here. It turns out I needed to use the Python library module cgi, and importantly, the text area needed an attribute: name='mytextarea'.

We're coding like it's 1995 :)

script.py

#!/usr/bin/python
import sys, os, cgi

s = '''Content-type: text/html

<head></head>'''

form = cgi.FieldStorage()

print s
print form.keys(),'<br>'
print form.getvalue('sex'),'<br>'
print form.getvalue('vehicle'),'<br>'
data = form.getvalue('mytextarea'),'<br>'
print len(data),'<br>'
print data[0]
print '</body></html>'

form.html

<Content-type: text/html>
<form name="input" action="cgi-bin/script.py" method="post">
First name: <input type="text" name="firstname" /><br />
Last name: <input type="text" name="lastname" /><br /><br />
<input type="radio" name="sex" value="male" /> Male<br />
<input type="radio" name="sex" value="female" /> Female<br /><br />
<input type="checkbox" name="vehicle" value="Bike" /> I have a bike<br />
<input type="checkbox" name="vehicle" value="Car" /> I have a car<br /><br />
<textarea name='mytextarea' rows="5" cols="40">type here</textarea>
<input type=submit value="Submit" />
</form>
</html>

Simple server running Python script on OS X: part 2

This continues the example from last time, using the built-in server on OS X and trying to understand how to transmit data from a form to a Python script. The form and the script files that I edit are on my Desktop, and after editing I have to remember to do:

cp script.py /Library/WebServer/CGI-Executables/script.py
cp form.html /Library/WebServer/Documents/form.html

I have Web Sharing turned on, with Apache set to only listen to localhost. I point Safari to:

http://localhost/form.html

and I get what's in the graphic. The form text is:

<Content-type: text/html>

<form name="input" action="cgi-bin/script.py" method="get">
    
First name: <input type="text" name="firstname" /><br />
Last name: <input type="text" name="lastname" /><br /><br />
<input type="radio" name="sex" value="male" /> Male<br />
<input type="radio" name="sex" value="female" /> Female<br /><br />
<input type="checkbox" name="vehicle" value="Bike" /> I have a bike<br />
<input type="checkbox" name="vehicle" value="Car" /> I have a car<br /><br />

<textarea rows="5" cols="40">
The cat was playing in the garden.
</textarea>
<input type="submit" value="Submit" />
</form>
</html>

Actually, in the source for this post I changed all the less-than symbols to < as usual. There has to be a blank line after the first one with Content-type, or I get an error.

The script is very simple. It uses os.environ to get values from the "environment":

#!/usr/bin/python
import os

s = '''Content-type: text/html

<head></head>'''

print s
D = os.environ
for k in D:
    print '<p>', k, D[k], '</p>'
print '</body></html>'

If I type data into the form and then hit submit, the form data goes to the script, which reads it and sends the result back as a web page:

HTTP_REFERER http://localhost/form.html

VERSIONER_PYTHON_PREFER_32_BIT no

SERVER_SOFTWARE Apache/2.2.15 (Unix) mod_ssl/2.2.15 OpenSSL/0.9.8l DAV/2

SCRIPT_NAME /cgi-bin/script.py

SERVER_SIGNATURE

REQUEST_METHOD GET

SERVER_PROTOCOL HTTP/1.1

QUERY_STRING firstname=Tom&lastname=Elliott&sex=male&vehicle=Bike&vehicle=Car

PATH /usr/bin:/bin:/usr/sbin:/sbin

HTTP_USER_AGENT Mozilla/5.0 (Macintosh; U; Intel Mac OS X 10_6_5; en-us) AppleWebKit/533.18.1 (KHTML, like Gecko) Version/5.0.2 Safari/533.18.5

HTTP_CONNECTION keep-alive

SERVER_NAME localhost

REMOTE_ADDR 127.0.0.1

SERVER_PORT 80

SERVER_ADDR 127.0.0.1

DOCUMENT_ROOT /Library/WebServer/Documents

SCRIPT_FILENAME /Library/WebServer/CGI-Executables/script.py

SERVER_ADMIN you@example.com

HTTP_HOST localhost

REQUEST_URI /cgi-bin/script.py?firstname=Tom&lastname=Elliott&sex=male&vehicle=Bike&vehicle=Car

HTTP_ACCEPT application/xml,application/xhtml+xml,text/html;q=0.9,text/plain;q=0.8,image/png,*/*;q=0.5

GATEWAY_INTERFACE CGI/1.1

REMOTE_PORT 51027

HTTP_ACCEPT_LANGUAGE en-us

__CF_USER_TEXT_ENCODING 0x46:0:0

VERSIONER_PYTHON_VERSION 2.6

HTTP_ACCEPT_ENCODING gzip, deflate

So, it looks pretty good.

There is a big thing missing, however. I do not yet know how to obtain the data from the textarea. Baby steps.

Simple server running Python script on OS X

This post is really simple. I went exploring the server that comes with stock OS X today, and these are my notes. What I'm thinking might be interesting is to use the browser as a poor man's GUI. I would have a page that lists available programs, that would link to a form to fill in any changes to the default values and the name of a file with the data. Then we'd run the analysis and write the file to disk, and launch Preview with the graphic.

What I accomplished today: running a Python script using cgi-bin that assembles a page which the server sends back, including data loaded from a file, and the result of a Python call to a datetime.date object.

At first I turned off Airport on my laptop before enabling Web Sharing in Sharing Prefs. That was too restrictive, so later I did as suggested here, I modified the config file:

sudo cp /etc/apache2/httpd.conf /etc/apache2/httpd.conf.orig
Password:
$ ls /etc/apache2
extra  httpd.conf.orig mime.types other
httpd.conf magic  original users

sudo cp /etc/apache2/httpd.conf ~/Desktop/x.txt

[UPDATE: Discussion on whether 127.0.0.1 can be "spoofed" and what good it would do here, here, here. Basically, no.]

But it wouldn't let me edit it at first! I didn't have the correct permissions.

localhost:Desktop telliott_admin$ ls -al x.txt
-rw-r--r--  1 root  staff  17727 Nov 16 13:11 x.txt
localhost:Desktop telliott_admin$ sudo chmod 661 x.txt
localhost:Desktop telliott_admin$ ls -al x.txt
-rw-rw---x  1 root  staff  17727 Nov 16 13:11 x.txt

The config file says:

# Change this to Listen on specific IP addresses as shown below to 
# prevent Apache from glomming onto all bound IP addresses.
#
#Listen 12.34.56.78:80
Listen 80

So I changed that line to Listen 127.0.0.1:80

sudo cp ~/Desktop/x.txt /etc/apache2/httpd.conf

I can point Safari at http://localhost/~telliott_admin/ or http://127.0.0.1/~telliott_admin/.
That's my "personal website" which serves this page: /Users/telliott_admin/Sites/index.html. The page has an image


<img src="images/gradient.jpg" alt="" 
height="304" width="800" border="0" />

The image is just where you'd expect it to be in /Users/telliott_admin/Sites/images

The page also has a couple of links


<a href="/manual/">Apache manual</a>
<a href="http://www.apache.org/httpd">Apache</a>

The manual referenced by the link is actually here: /Library/WebServer/share/httpd/manual/index.html

My "computers website" is: http://localhost/

This page just says "It works!" It is here: /Library/WebServer/Documents/index.html.en

So I guess when the URL ends with '/' you load the index.html page

There is also an empty directory /Library/WebServer/CGI-Executables, which we're going to use next.

I made a file called data.txt on my Desktop:

No point without some data.
This is my data.

cp data.txt /Library/WebServer/Documents/data.txt

I made a second file with a Python script, called script.py:

#!/usr/bin/python
import os,time
from datetime import date
today = date.today()

fn = '/Library/WebServer/Documents/data.txt'
#fn = '/Documents/data.txt'
FH = open(fn,'r')
data = FH.read()
FH.close()

s = '''Content-type: text/html

<head>Title</head>
<hr>
<p></p>
<body>xyz
<hr>
<p></p>'''

print s, data,
print '<br>'
print 'and this is my date  :)'
print '<br>'
print today
print '</body></html>'

$ cp script.py /Library/WebServer/CGI-Executables/script.py
$ sudo chmod +x /Library/WebServer/CGI-Executables/script.py
Password:
$ ls -al /Library/WebServer/CGI-Executables/script.py
-rwxrwxrwx@ 1 telliott_admin  staff  298 Nov 16 13:37 /Library/WebServer/CGI-Executables/script.py

Point browser at: http://localhost/cgi-bin/script.py

When it doesn't work, do this:

$ cat /var/log/apache2/error_log

Monday, November 15, 2010

Lucretius speaks to us of atoms

Clothes hung above a surf-swept shore grow damp; spread in the sun they dry again. Yet it is not apparent to us how the moisture clings to the cloth, or flees the heat. Water, then, is dispersed in particles, atoms too small to be observable.
..
For surely the atoms did not hold council assigning order to each, flexing their keen minds with questions of place and motion and who goes where. But shuffled and jumbled in many ways, in the course of endless time they are buffeted, driven along, chancing upon all motions, combinations. At last they fall into such an arrangement as would create this universe..

-Lucretius (De rerum natura)
as quoted by David Lindley in Boltzmann's Atom

Boltzmann by Lindley

A very nice book about the history of atoms, and Boltzmann. Although it was out of print, I bought it through abebooks for $$$, because I'd read Lindley's other books. This one is great. And the notation on the withdrawal stamp, simply sad ("out-of-date, no longer needed" ??). It is virtually unmarked.

Installing SciPy with MacPorts

To continue with the previous post (here), I thought it would make life easier if I could install SciPy or other software using MacPorts. Some years have passed since the unpleasantness with Fink (so long that I'm a little hazy now on the details). Anyway, as explained previously, after installing MacPorts, I ran into a problem with not being able to build one of SciPy's dependencies: LAPACK.

MacPorts uses a file called a Portfile to control all of what happens when a package is downloaded and built. The problem stems from the release of a new version of LAPACK. Apparently they didn't tell MacPorts about it, so the Portfile for the atlas port (which includes LAPACK), is out of date. In particular, the Portfile directions allow verification of three checksums on the download (md5, sha-1 and rmd160). Since the Portfile was out of date these didn't match. You can do this manually on any file of interest like so:

md5 filename.tar.gz
openssl sha1 filename.tar.gz
openssl rmd160 filename.tar.gz

The way to fix this issue is to change the Portfile. For example, the MacPorts stuff is mostly or all in

$ ls /opt/local/var/macports
build  logs  port-help.tcl registry sources
distfiles packages receipts software

and the Portfile is there:

ls /opt/local/var/macports/sources/rsync.macports.org/release/ports/math/atlas
Portfile Portfile.orig Portfile.rej files

A standard approach would be to verify that the download obtained by MacPorts is authentic, by verifying the checksums against a download you obtain yourself, or values published by a trusted source, then edit the MacPorts file to change the values we're checking against. One could simply copy this modified Portfile and replace the original. However, the correct method is a little different, allowing one to easily reverse the patch if necessary later.

The "manual" method creates a Portfile Patch by doing a diff between the copy and the original, and then calling patch:

cp -p Portfile Portfile.orig
diff -u Portfile.orig Portfile > Portfile-atlas-fixed.diff
patch -p0 < ~/Desktop/Portfile-atlas-fixed.diff

However, then I ran into a second problem, which was that the new version of LAPACK was extracted and configured but there was problem later in the process (it wasn't actually in the build phase but (perhaps) in clean). Since the LAPACK release is so new, I guess that nothing will depend on it being version 3.3.0, though a lot has changed. Rather than figure out what the issue is, I edited the Portfile in a different way, keeping the old checksums, and changing the download we request. That's the screenshot at the top of the post. Textmate knows a diff file when it see one, and syntax colors appropriately. That answers the question I posed at the end last time (here).

By the way, the basic reason for the checksum mismatch is that the standard LAPACK download file does not have a version number, at least the tgz doesn't. If we'd asked for the download by version number, we'd have got the correct file. The "stealth upgrade" is because the default download is unnumbered.

Whether I did the direct copy or the patch thing on the version that worked, I don't exactly remember now. patch has an issue that it wants to know whether a previous patch needs to be reversed or not, and can fail. Not sure now, but a direct copy will certainly work.

Now we use a nifty trick involving a shell variable (I think that's what it's called) to cd into where we need to be, and finish what we're trying right now:

$(port dir atlas)
-bash: /opt/local/var/macports/sources/rsync.macports.org/release/ports/math/atlas: is a directory

cd $(port dir atlas)
sudo cp ~/Desktop/Portfile Portfile
c-98-236-78-154:atlas telliott_admin$ sudo port install atlas
Portfile changed since last build; discarding previous state.
--->  Computing dependencies for atlas
--->  Fetching atlas
--->  Verifying checksum(s) for atlas
--->  Extracting atlas
--->  Applying patches to atlas
--->  Configuring atlas
--->  Building atlas
--->  Staging atlas into destroot
--->  Installing atlas @3.8.3_4+gcc44
--->  Activating atlas @3.8.3_4+gcc44
--->  Cleaning atlas
c-98-236-78-154:atlas telliott_admin$

Yes!

sudo port install py26-scipy @0.8.0
..
--->  Fetching py26-scipy
--->  Attempting to fetch scipy-0.8.0.tar.gz from http://superb-east.dl.sourceforge.net/scipy
--->  Attempting to fetch scipy-0.8.0.tar.gz from http://downloads.sourceforge.net/scipy
--->  Verifying checksum(s) for py26-scipy
--->  Extracting py26-scipy
--->  Configuring py26-scipy
--->  Building py26-scipy
--->  Staging py26-scipy into destroot
--->  Installing py26-scipy @0.8.0_0+gcc44
--->  Activating py26-scipy @0.8.0_0+gcc44
--->  Cleaning py26-scipy

$ /opt/local/bin/python2.6
Python 2.6.6 (r266:84292, Nov 14 2010, 17:37:27) 
[GCC 4.2.1 (Apple Inc. build 5664)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from scipy import *
>>>

How cool is that?

Sunday, November 14, 2010

MacPorts

I decided to try MacPorts, but ran into some problems that aren't solved yet. I thought maybe it would make it easier to install SciPy. (Yes, I know I've already got it installed). I used the MacPorts installer first and then tried working on Python first and afterwards SciPy.

MacPorts wants you to keep everything separate, so for example, I need a different Python installed by them. The first weird thing is they came up with numerous dependencies for Python, which I never heard of. One of those (db46) failed because I didn't have the very latest Java installed. So I (re) registered as a Developer, and got that from Apple:

javadeveloper_10.6_10m3261.dmg

Next:

sudo port install python27
Password:
--->  Computing dependencies for python27
--->  Dependencies to be installed: db46 gdbm gettext expat libiconv gperf ncurses ncursesw openssl zlib readline sqlite3 tk Xft2 fontconfig freetype pkgconfig xrender xorg-libX11 xorg-bigreqsproto xorg-inputproto xorg-kbproto xorg-libXau xorg-xproto xorg-libXdmcp xorg-libxcb python26 xorg-libpthread-stubs xorg-xcb-proto libxml2 xorg-util-macros xorg-xcmiscproto xorg-xextproto xorg-xf86bigfontproto xorg-xtrans xorg-renderproto tcl xorg-libXScrnSaver xorg-libXext xorg-scrnsaverproto

Wow. Scroll that out there to see them all. It's a long list.

Two strange things: in the middle of doing Python 2.7 it downloaded and built Python 2.6! And...

--->  Fetching gcc44
--->  Attempting to fetch gcc-core-4.4.5.tar.bz2 from http://mirror.facebook.net/gnu/gnu/gcc/gcc-4.4.5
--->  Attempting to fetch gcc-fortran-4.4.5.tar.bz2 from http://mirror.facebook.net/gnu/gnu/gcc/gcc-4.4.5
--->  Attempting to fetch gcc-g++-4.4.5.tar.bz2 from http://mirror.facebook.net/gnu/gnu/gcc/gcc-4.4.5
--->  Attempting to fetch gcc-java-4.4.5.tar.bz2 from http://mirror.facebook.net/gnu/gnu/gcc/gcc-4.4.5
--->  Attempting to fetch gcc-objc-4.4.5.tar.bz2 from http://mirror.facebook.net/gnu/gnu/gcc/gcc-4.4.5
--->  Verifying checksum(s) for gcc44
--->  Extracting gcc44
--->  Configuring gcc44
--->  Building gcc44
--->  Staging gcc44 into destroot
--->  Installing gcc44 @4.4.5_0
--->  Activating gcc44 @4.4.5_0
--->  Cleaning gcc44

(It took a very long time). What is harder to believe? That Facebook is a mirror for gnu, or that these guys want to install GCC 4.4.5? Well, they did. Next up: SciPy. And of course, there is no package for Python 2.7, only 2.6!

sudo port install py26-scipy @0.8.0
Password:
--->  Computing dependencies for py26-scipy
--->  Dependencies to be installed: atlas py26-nose py26-distribute py26-numpy fftw-3 swig-python bison m4 gsed python_select swig pcre
--->  Verifying checksum(s) for atlas
Error: Checksum (md5) mismatch for lapack.tgz
Error: Checksum (sha1) mismatch for lapack.tgz
Error: Checksum (rmd160) mismatch for lapack.tgz
Error: Target org.macports.checksum returned: Unable to verify file checksums
Error: Failed to install atlas

A really stupid but quite common problem. The md5 checksum in the portfile equals the md5 checksum of the previous version of lapack(3.2.2). So I started reading the guide, and figured out how to deal with the issue. I downloaded the original software using curl to maintain it as a .tgz:

curl -O http://www.netlib.org/lapack/lapack.tgz

computed various checksums

md5 lapack.tgz 
MD5 (lapack.tgz) = 70aba8e5ecdccb6003850db178e551a2

openssl sha1 lapack.tgz 
SHA1(lapack.tgz)= a0354c8eda9737319f93472068bbf187b26e1e69

openssl rmd160 lapack.tgz 
RIPEMD160(lapack.tgz)= 253d0597f275fd5cd86a1a447ef56c92635aff74

Followed the instructions:

Make a copy of the Portfile you wish to modify; both files must be in the same directory, though it may be any directory.

cd $(port dir atlas)
cp Portfile ~/Desktop
cp Portfile ~/Desktop/Portfile.orig

cd ~/Desktop
diff -u Portfile.orig Portfile > Portfile-atlas.diff
cd $(port dir atlas)
sudo patch -p0 < ~/Desktop/Portfile-atlas.diff
Password:
patching file Portfile

ls -al Portfile 
-rw-r--r--  1 root  wheel  11610 Nov 15 06:04 Portfile

And now MacPorts unpacks the thing, but the build fails later:

sudo port install atlas
Portfile changed since last build; discarding previous state.
--->  Computing dependencies for atlas
--->  Fetching atlas
--->  Verifying checksum(s) for atlas
--->  Extracting atlas
--->  Applying patches to atlas
--->  Configuring atlas
--->  Building atlas
Error: Target org.macports.build returned: shell command failed (see log for details)
Log for atlas is at: /opt/local/var/macports/logs/_opt_local_var_macports_sources_rsync.macports.org_release_ports_math_atlas/main.log
Error: Status 1 encountered during processing.
To report a bug, see <http://guide.macports.org/#project.tickets>

I haven't solved that yet. And that seems to be a weakness of what looks like a pretty sophisticated system. How would you specify that it use a previous version of LAPACK?

New plotter for phylogenetic trees: version 0.1

This is the last post on this project. The previous posts should be listed in the archive in the sidebar.

I have a phylogenetics project for which I constructed a boutique database. The definition is in the file db.groups.txt, one section of which looks like this:


aeromonas
X60415.1       Aeromonas_trota_ATCC49657
X74677.1       Aeromonas_hydrophila_ATCC7966T
EU770300.1     Aeromonas_enteropelogenes_MS12

set off from surrounding entries by double newlines. The "group" aeromonas has three sequences. There are also super-groups that correspond to bacterial Phyla or sub-Phyla, e.g., the gamma-Proteobacteria. These are defined in the file big.groups.txt like this:


gamma
aeromonas
cardio
pseudo
moraxella
entero
haemo1
haemo2
steno
xantho

Normally the sequences would be fetched from Genbank by a slightly complicated script, but I'd like you to be able to follow along if you want, so we'll do it almost manually. The file fetchSeqs.py contains a function that parses these two files and grabs the info we need for the "gamma" super-group. It also does a manual fetch from Genbank. We save the data to disk like so:

python fetchSeqs.py > seqs.txt

The beginning of one of the entries looks like this:

>gi|312967|emb|X60415.1| Aeromonas trota 16S rRNA gene, strain ATCC 49657
GAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAGCGGCAGCGGGAA

As you can see, the title line of this FASTA-formatted text is really long. So in the second script, analyzeSeqs.py, we use the same info from before to replace the long title with something shorter. This script also uses MUSCLE to align the sequences. You need to have it installed for this to work. Or you could use Clustal. I used the ape library in R to make a neighbor-joining tree and write this to disk as 'tree.txt'. I show a plot of what R gave me (actually the real process for this one was more complicated, resulting the pretty colors). But at least you can see what the tree is supposed to look like.

The last module is test_plotter.py It shows how to make a couple different kinds of plots with our plotter. All of them have colored node labels as defined in the script. The first plot is the standard one, the graphic is at the top of the post. In the next, we plot with the internal node labels showing, so that we can identify the name of a node to re-root the tree, if wanted.

In the third example, I got a balanced tree from PyCogent and plot that.

You might notice that the font used for the external node labels is now italic as it should be. It doesn't look that hot, I think that matplotlib is not using the OS X system fonts. Not sure why that is, yet.

I put a zip of all the files up on Dropbox (here).

New plotter for phylogenetic trees: customization

Continuing with a tree-plotter in Python. Previous posts here, here, and here. These examples depend on PyCogent and matplotlib.

I've got a somewhat bigger tree here, using sequences from a project we're working on. I used clustal and R to make a neighbor-joining tree. As you can see in the code, there is a lot of flexibility in terms of customized colors, etc.

code listing:

import tree_utils as TU
import tree_plotter as TP

fn = 'veillonella.tree.txt'
D = TU.make_tree_dict(fn=fn)
tr = D['meta']['tree']

attr = TP.get_default_attr()
attr['e_node_dot_size'] = 35
attr['figure_type'] = 'png'
attr['e_node_label_size'] = 10
attr['e_node_bar_color'] = 'k'

attr['using_alternate_names'] = True
e_node_names = D['meta']['e_node_names']
#L = [e.replace('_',' ') for e in e_node_names]
L = [e.split('_')[-1] for e in e_node_names]
# lose the version numbers on Genbank seqs
L = [e.split('.')[0] for e in e_node_names]
alt_name_dict = dict(zip(e_node_names,L))
attr['alternate_name_dict'] = alt_name_dict

attr['using_e_node_specific_colors'] = True
def f(name):
    if 'DB' in name:  return 'b'
    return 'r'
L = [f(e) for e in e_node_names]
node_label_color_dict = dict(zip(e_node_names,L))
attr['node_label_color_dict'] = node_label_color_dict

def f(name):
    if 'DB' in name:  return 'b'
    return 'r'
L = [f(e) for e in e_node_names]
dot_color_dict = dict(zip(e_node_names,L))
attr['dot_color_dict'] = dot_color_dict

TP.plot(D,attr=attr)

'''
Rcode:
library(ape)
setwd('Desktop')
dna = read.dna('veillonella.fasta',format='fasta')
tree = nj(dist.dna(dna))
plot(tree,cex=1)
axisPhylo()
write.tree(tree,'veillonella.tree.txt')
'''

New plotter for phylogenetic trees: re-rooting

Continuing with a tree-plotter in Python. Previous posts here and here. These examples depend on PyCogent and matplotlib.

One of the great things about building on PyCogent is we can use its re-rooting routine. In the output for this example we print the same tree as before, rooted at all possible internal nodes.

I've made a few tweaks to the code. make_tree_dict will take either a tree_string, a filename or a PyCogent tree. For the graphic above, we grab the default plot attributes, make the internal nodes visible, and then call plot. The code for this example is at the bottom. (I'm not going to post modifications to the other code until it's finalized. If you want it, shoot me an email).

To do list:
• learn matplotlib typography (e.g. italics)
• adjust for label sizes in a smart way
• provide for printing bootstrap values or internal node names

output:

$ python test_plotter.py 
original root
                    /-Stenotrophomonas_maltophilia
          /edge.0--|
         |          \-Kingella_oralis
         |
-root----|--Pseudomonas_aeruginosa
         |
         |                    /-Salmonella_typhi
         |          /edge.1--|
          \edge.2--|          \-Escherichia_coli
                   |
                    \-Haemophilus_parainfluenzae
------------------------------------------------------------
edge.0
          /-Stenotrophomonas_maltophilia
         |
         |--Kingella_oralis
-root----|
         |          /-Pseudomonas_aeruginosa
         |         |
          \edge.0--|                    /-Salmonella_typhi
                   |          /edge.1--|
                    \edge.2--|          \-Escherichia_coli
                             |
                              \-Haemophilus_parainfluenzae
------------------------------------------------------------
edge.2
                    /-Salmonella_typhi
          /edge.1--|
         |          \-Escherichia_coli
         |
-root----|--Haemophilus_parainfluenzae
         |
         |                    /-Stenotrophomonas_maltophilia
         |          /edge.0--|
          \edge.2--|          \-Kingella_oralis
                   |
                    \-Pseudomonas_aeruginosa
------------------------------------------------------------
edge.1
          /-Salmonella_typhi
         |
         |--Escherichia_coli
-root----|
         |          /-Haemophilus_parainfluenzae
         |         |
          \edge.1--|                    /-Stenotrophomonas_maltophilia
                   |          /edge.0--|
                    \edge.2--|          \-Kingella_oralis
                             |
                              \-Pseudomonas_aeruginosa
------------------------------------------------------------
balanced tree
                    /-Stenotrophomonas_maltophilia
          /edge.0--|
         |          \-Kingella_oralis
         |
-root----|--Pseudomonas_aeruginosa
         |
         |                    /-Salmonella_typhi
         |          /edge.1--|
          \edge.2--|          \-Escherichia_coli
                   |
                    \-Haemophilus_parainfluenzae

code listing:

import tree_utils as TU
import tree_plotter as TP

fn = 'tree.txt'
D = TU.make_tree_dict(fn=fn)
tr = D['meta']['tree']
print 'original root'
print tr.asciiArt()
print '-'*60

L = D['meta']['i_node_names']
for name in L[1:]:
    print name
    print tr.rootedAt(name).asciiArt()
    D = TU.make_tree_dict(tr=tr)
    print '-'*60

print 'balanced tree'
print tr.balanced().asciiArt()

tr = tr.rootedAt('edge.2')
D = TU.make_tree_dict(tr=tr)
attr = TP.get_default_attr()
attr['i_node_dots_visible'] = True
TP.plot(D,attr=attr,ofn='balanced')

Saturday, November 13, 2010

New plotter for phylogenetic trees: plotting

Here is a first pass at a plotter, following up on the last post.

~~It's ugly code yet, but~~ I thought I would show you what I have. What I need to do is really think harder about how this module should be organized. But you can compare the result with what we got here.

[UPDATE: Modified the code to be easily customizable. Next up: more extensive testing.]
UPDATE2: Modified the code yet again.]

output:

         alternate name dict   None
    default e node dot color   k
              dot color dict   None
            e node bar color   r
       e node bar color list   None
             e node dot size   75
         e node dots visible   True
  e node label default color   k
           e node label size   14
       e node labels visible   True
                 figure type   png
     horizontal axis visible   True
            i node bar color   k
            i node dot color   magenta
         i node dots visible   False
   i node vertical bar color   k
          label width factor   1.1
                  line width   2
       node label color dict   None
            r node dot color   orange
       using alternate names   False
using e node specific colors   False
       vertical axis visible   False

code listing:

import matplotlib.pyplot as plt
import tree_utils as tu

def plot(D,attr=None,ofn=None):
    if not attr:
        attr = get_default_attr()
    SZ = attr['e_node_dot_size']
    lw = attr['line_width']
    maxx, maxy = D['meta']['max_xy']
    
    all_node_names = D['meta']['all_node_names']
    e_node_names = D['meta']['e_node_names']
    i_node_names = D['meta']['i_node_names']

    # actual node dictionaries
    e_node_dicts = [D[n] for n in e_node_names]
    i_node_dicts = [D[n] for n in i_node_names]

    # extract the values we need
    e_node_positions = [(nD['x'],nD['y']) for nD in e_node_dicts]
    e_node_lengths = [(nD['dist_to_parent']) for nD in e_node_dicts]
    i_node_positions = [(nD['x'],nD['y']) for nD in i_node_dicts]
    i_node_lengths = [(nD['dist_to_parent']) for nD in i_node_dicts]
    i_node_verticals = [(nD['y_bott'], nD['y_top']) for nD in i_node_dicts]

    # external nodes
    for i in range(len(e_node_names)):
        x,y = e_node_positions[i]
        d = e_node_lengths[i]
        xleft = x - d
        # bars
        L = attr['e_node_bar_color_list']
        if L:
            c = L[i]
        else:
            c = attr['e_node_bar_color']
        plt.plot([xleft,x],[y,y],color=c,lw=lw,zorder=1)
        # dots
        name = e_node_names[i]
        if attr['e_node_dots_visible']:
            if attr['using_e_node_specific_colors']:
                c = attr['dot_color_dict'][name]
            else:
                c = attr['default_e_node_dot_color']
        plt.scatter(x,y,color=c,s=SZ,zorder=2)
    
    # internal nodes, root is i = 0
    for i in range(len(i_node_names))[1:]:
        x,y = i_node_positions[i]
        d = i_node_lengths[i]
        y_bott, y_top = i_node_verticals[i]
        x0 = x - d
        # bars
        c = attr['i_node_bar_color']
        plt.plot([x0,x],[y,y],color=c,lw=lw,zorder=1)
        # verticals
        c = attr['i_node_vertical_bar_color']
        plt.plot([x,x],[y_bott, y_top],color=c,lw=lw,zorder=1)
        # dots
        if attr['i_node_dots_visible']:
            c = attr['i_node_dot_color']
            plt.scatter(x,y,color=c,s=SZ,zorder=2)
    
    # root
    i = 0
    x,y = i_node_positions[i]
    y_bott, y_top = i_node_verticals[i]
    # verticals
    c = attr['i_node_vertical_bar_color']
    plt.plot([x,x],[y_bott, y_top],color=c,lw=lw,zorder=1)
    # dots
    if attr['i_node_dots_visible']:
        c = attr['r_node_dot_color']
        plt.scatter(x,y,color=c,s=SZ,zorder=2)

    L = [len(name) for name in e_node_names]
    max_label_width = maxx/100.0 * max(L)
    max_label_width *= attr['label_width_factor']

    # external node labels
    if attr['e_node_labels_visible']:
        for i in range(len(e_node_names)):
            x,y = e_node_positions[i]
            dx = maxx/100.0 * 4
            name = e_node_names[i]
            if attr['using_alternate_names']:
                s = attr['alternate_name_dict'][name]
            else:
                s = name
            if attr['using_e_node_specific_colors']:
                c = attr['node_label_color_dict'][name]
            else:
                c = attr['e_node_label_default_color']
            plt.text(x + dx, y, s,
                fontname = 'Helvetica', 
                fontsize = attr['e_node_label_size'],
                color = c,
                ha = 'left',va = 'center')

    ax = plt.axes()
    if not attr['vertical_axis_visible']:
        ax.yaxis.set_visible(False)
    if not attr['horizontal_axis_visible']:
        ax.xaxis.set_visible(False)
    ax.set_xlim(-maxx/10.0,maxx*1.1 + max_label_width)
    ax.set_ylim(-maxy/10.0,maxy*1.1)
    
    if ofn:  
        plt.savefig(ofn + '.' + attr['figure_type'])
    else:
        plt.savefig('example.' + attr['figure_type'])
    
def get_default_attr():
    attr = dict()
    attr['e_node_dot_size'] = 75
    attr['line_width'] = 2
    attr['figure_type'] = 'png'
    attr['horizontal_axis_visible'] = True
    attr['vertical_axis_visible'] = False
    
    attr['e_node_bar_color'] = 'r'
    attr['e_node_bar_color_list'] = None
    attr['e_node_dots_visible'] = True
    attr['e_node_labels_visible'] = True
    attr['default_e_node_dot_color'] = 'k'
    attr['e_node_label_size'] = 14
    attr['e_node_label_default_color'] = 'k'
    attr['label_width_factor'] = 1.1
    
    attr['i_node_bar_color'] = 'k'
    attr['i_node_vertical_bar_color'] = 'k'
    attr['i_node_dots_visible'] = False
    attr['i_node_dot_color'] = 'magenta'
    attr['r_node_dot_color'] = 'orange'
    
    attr['using_alternate_names'] = False
    attr['alternate_name_dict'] = None
    attr['using_e_node_specific_colors'] = False
    attr['node_label_color_dict'] = None
    attr['dot_color_dict'] = None
    return attr
    
def print_defaults():
    attr = get_default_attr()
    N = max([len(k) for k in attr.keys()])
    for k in sorted(attr.keys()):
        v = attr[k]
        k = k.replace('_',' ')
        print k.rjust(N), ' ', v

if __name__ == '__main__':
    fn = 'tree.txt'
    file_data = tu.load_data(fn)
    D = tu.make_tree_dict(ts=file_data)
    plot(D)
    print_defaults()

New plotter for phylogenetic trees: parsing

In July of last year, I had a series of posts about drawing phylogenetic trees (here, here, here, here and here). Some very early stuff is here and here.

I want to revisit this topic to try to get a tree plotter in python that is really flexible. The important thing to remember here is that we do not want to "roll our own" phylogenetics software. We should let the experts do that for us (and debug it!). Now that we know about PyCogent (web page, links in the sidebar), I say let them do the heavy lifting.

So I wrote a module today that starts with a Newick format phylogenetic tree (the one we used before here), and gets PyCogent to figure out all its connections. You can see how versatile the PyCogent code for dealing with trees is in the docs.

This is the tree:

((Stenotrophomonas_maltophilia:0.07574,
Kingella_oralis:0.08026)
:0.00827,
Pseudomonas_aeruginosa:0.05950,
((Salmonella_typhi:0.01297,
Escherichia_coli:0.01491)
:0.03356,
Haemophilus_parainfluenzae:0.06113)
:0.03863);

Anyway, here is the output from the module. It's not very well tested yet. So I hope you find a bug and let me know. But I'm not Donald Knuth so I'm not going to pay you :)

[UPDATE: I replaced the code with a later version]

output:

$ python tree_utils.py
making default tree 

finding children
edge.0
['Stenotrophomonas_maltophilia', 'Kingella_oralis']
edge.1
['Salmonella_typhi', 'Escherichia_coli']
edge.2
['edge.1', 'Haemophilus_parainfluenzae']
root
['edge.0', 'Pseudomonas_aeruginosa', 'edge.2']

Stenotrophomonas_maltophilia
node 'Stenotrophomonas_maltophilia':0.07574;
ancestors ['edge.0', 'root']
name Stenotrophomonas_maltophilia
is_root False
y 0.0
x 0.08401
dist_to_parent 0.07574
is_external True

Kingella_oralis
node 'Kingella_oralis':0.08026;
ancestors ['edge.0', 'root']
name Kingella_oralis
is_root False
y 0.019952
x 0.08853
dist_to_parent 0.08026
is_external True

Pseudomonas_aeruginosa
node 'Pseudomonas_aeruginosa':0.0595;
ancestors ['root']
name Pseudomonas_aeruginosa
is_root False
y 0.039904
x 0.0595
dist_to_parent 0.0595
is_external True

Salmonella_typhi
node 'Salmonella_typhi':0.01297;
ancestors ['edge.1', 'edge.2', 'root']
name Salmonella_typhi
is_root False
y 0.059856
x 0.08516
dist_to_parent 0.01297
is_external True

Escherichia_coli
node 'Escherichia_coli':0.01491;
ancestors ['edge.1', 'edge.2', 'root']
name Escherichia_coli
is_root False
y 0.079808
x 0.0871
dist_to_parent 0.01491
is_external True

Haemophilus_parainfluenzae
node 'Haemophilus_parainfluenzae':0.06113;
ancestors ['edge.2', 'root']
name Haemophilus_parainfluenzae
is_root False
y 0.09976
x 0.09976
dist_to_parent 0.06113
is_external True

root
y_top 0.084796
ancestors []
name root
is_root True
immediate_children ['edge.0', 'Pseudomonas_aeruginosa', 'edge.2']
y_bott 0.009976
y 0.039904
x 0
dist_to_parent 0
is_external False

edge.0
y_top 0.019952
ancestors ['root']
name edge.0
is_root False
immediate_children ['Stenotrophomonas_maltophilia', 'Kingella_oralis']
y_bott 0.0
y 0.009976
x 0.00827
dist_to_parent 0.00827
is_external False

edge.2
y_top 0.09976
ancestors ['root']
name edge.2
is_root False
immediate_children ['edge.1', 'Haemophilus_parainfluenzae']
y_bott 0.069832
y 0.084796
x 0.03863
dist_to_parent 0.03863
is_external False

edge.1
y_top 0.079808
ancestors ['edge.2', 'root']
name edge.1
is_root False
immediate_children ['Salmonella_typhi', 'Escherichia_coli']
y_bott 0.059856
y 0.069832
x 0.07219
dist_to_parent 0.03356
is_external False

code listing for tree_utils.py:

import sys
from cogent import LoadTree

s = '''
((Stenotrophomonas_maltophilia:0.07574,
Kingella_oralis:0.08026)
:0.00827,
Pseudomonas_aeruginosa:0.05950,
((Salmonella_typhi:0.01297,
Escherichia_coli:0.01491)
:0.03356,
Haemophilus_parainfluenzae:0.06113)
:0.03863);'''

def load_data(fn):
    FH = open(fn,'r')
    data = FH.read().strip()
    FH.close()
    return data

# for flexibility, we can this function with
# a tree_string
# a filename
# or a PyCogent tree
def make_tree_dict(ts=None,fn=None,tr=None,debug=False):
    if tr:
        pass
    elif fn:
        tr = LoadTree(filename=fn)
    else:
        tr = LoadTree(treestring=ts)
    nodes = tr.getNodesDict()
    
    # and their names
    all_node_names = tr.getNodeNames()
    # PyCogent organizes the external node names
    # in the correct order, so the rest is easy

    # rearrange the data my way
    D = dict()
    D['meta'] = { 'all_node_names':all_node_names }
    D['meta']['tree'] = tr
    e_node_names = list()
    i_node_names = list()
    
    for name in all_node_names:
        node = nodes[name]
        external = node.isTip()
        if external:  e_node_names.append(name)
        else:  i_node_names.append(name)
        
        # in order up the tree
        aL = [a.Name for a in node.ancestors()]
        dist_to_root = node.distance(nodes['root'])
        
        nD = {'name':name,
              'node':node,
              'is_external':external,
              'is_root': name == 'root',
              'x':dist_to_root,
              'ancestors':aL }
              
        if not nD['is_root']:
            dist = node.distance(nodes[aL[0]])
            nD['dist_to_parent'] = dist
        else:
            nD['dist_to_parent'] = 0
        D[name] = nD
    
    D['meta']['e_node_names'] = e_node_names
    D['meta']['i_node_names'] = i_node_names
    D = compute_y_pos(D,debug=debug)
    if debug:  show(D)
    return D
    
def compute_y_pos(D,debug=False):
    all_node_names = D['meta']['all_node_names']
    e_node_names = D['meta']['e_node_names']
    i_node_names = D['meta']['i_node_names']

    # actual node dictionaries
    e_node_dicts = [D[n] for n in e_node_names]
    i_node_dicts = [D[n] for n in i_node_names]

    # just make a square plot
    xmax = max([nD['x'] for nD in e_node_dicts])
    ymax = xmax
    D['meta']['max_xy'] = xmax, ymax

    # compute y for external nodes
    N = 1.0*(len(e_node_names)-1)
    for nD in e_node_dicts:
        i = e_node_names.index(nD['name'])
        f = i/N
        nD['y'] = f*ymax

    def mean(L):  return sum(L)*1.0/len(L)

    def get_child_node_names(parent):
        if debug:  print parent
        cL = list()
        for name in all_node_names:
            if name == 'root':  continue
            nD = D[name]
            if parent == nD['ancestors'][0]:
                cL.append(name)
        if debug:  print cL
        return cL

    # compute y for internal nodes
    # must fill these from the bottom up
    L = i_node_names[:]
    #L.sort()     # sorts 0,1,2 etc. then root
    def f(name):  # sort by x-value !
        nD = D[name]
        return nD['x']
    L = sorted(L,key=f,reverse=True)
    
    if debug:  print 'finding children'    
    for name in L:
        nD = D[name]
        cL = get_child_node_names(parent=name)
        nD['immediate_children'] = cL
        yL = [D[child]['y'] for child in cL]
        y0 = min(yL)
        y1 = max(yL)
        if len(yL) == 3:
            y = sorted(yL)[1]
        else:
            y = mean([y0,y1])        
        nD = D[name]
        nD['y'] = y
        nD['y_top'] = y1
        nD['y_bott'] = y0
    if debug:  print
    
    # sanity check for ancestor lists
    for name in all_node_names:
        if name == 'root':  continue
        assert D[name]['ancestors'][-1] == 'root'
    return D
    
def show(D):
    e_node_names = D['meta']['e_node_names']
    i_node_names = D['meta']['i_node_names']
    for name in e_node_names + i_node_names:
        print name
        nD = D[name]
        for k in nD:
            if k == 'node' and not nD['is_external']: 
                continue
            print k, nD[k]
        print

if __name__ == '__main__':
    tree_string = None
    if len(sys.argv) > 1:
        fn = sys.argv[1]
        try:
            tree_string = load_data(fn)
        except IOError:
            print 'could not open file:', fn
            sys.exit()
    if not tree_string:
        print 'making default tree', '\n'
        tree_string = s.strip()
    D = make_tree_dict(ts=tree_string,debug=True)

Moses wrote the U.S. Constitution

wikimedia

Read about it here. All I can say is, holy sh*t!

Things to explore

This looks very interesting:

As does this:

Correction term for the mean

A few days ago I asked a question on Stack Exchange about the "correction term for the mean" and two different ways of calculating the sum of squares for variance (original post here).

The answer as formulated by Srikant Vadali is developed elegantly in the first link, but recapping with respect to the math, the answer is pretty easy when you know how ;-)

If X is an array of numbers and

n = len(X)
m = sum(X)/n

written another way:

m n = sum(X)

then expanding the original form of the sum of squares:

sum(X - m)² = sum( X² - 2 m X + m² )
       = sum(X²) - sum(2 m X) + sum(m²)

But m is a constant so it can move in front of the summation. The second term is then

-sum(2 m X) = -2 m sum(X)
    = -2 m² n

And the third term is

sum(m²) = m² n

They add to give

- m² n

but now we can go back to sum(X):

- m sum(X)

and the whole thing is

sum(X - m)² = sum(X²) - m sum(X)

which was "precisely what was required to be proved." (here)

I think the existence of the family of sites including Stack Exchange and Stack Overflow is just fantastic.

More interesting than my transcription of Srikant's derivation is that his markup looks so much prettier. The source for the first line is

$\sum_i(X_i-m)^2 = \sum_i(X_i^2 + m^2 - 2 X_i m)$

which is obviously not rendering correctly on this page but looks a bit like LaTeX. Anybody know how this is done or how one could do it in Blogger? Time for another question, I guess.

Silent sites in EMBOSS

EMBOSS includes a program called silent that determines positions within a coding sequence which can be mutated to give a convenient restriction site without changing the encoded protein.

$ silent hemA.txt
Find restriction sites to insert (mutate) with no translation change
Comma separated enzyme list [all]: 
Output report [hema.silent]: emboss.results.txt

The column headings and a sample line are as shown, when run on the sequence of the Salmonella typhimurium hemA gene:

  Start     End  Strand EnzymeName RS-Pattern     Base-Posn    AAs Silent Mutation
    327     332       + SpeI       ACTAGT               330    L.L    Yes     G->A

In an early post I blogged about this and posted some code.

It seemed like a good idea to compare results from the two programs. One problem is the redundancy of restriction enzymes with respect to cleavage sites. Every position reported by silent in EMBOSS has a number of enzymes, some as many two dozen or so (each one listed as a separate hit).

If you've ever worked with these enzymes you know that they have personalities, some are highly active, some produce overhangs, some are stable, some not so good. Filtering out the redundancy in a smart way is a bit of a pain. That's why one of the optional inputs to the program is a list of enzymes that you like. I ran silent with the default setting (lots of hits), and then for this analysis I made a short list of "good_enzymes."

To repeat what we did previously I copied out Restriction_Dictinary.py from the biopython-1.55 source. Then I went back to the old post and got three files (REnzymes, GeneticCode2, extrasites).

I ran REnzymes.py and it looks like it works even though the format of Restriction_Dictinary.py has changed (and is really unrecognizable to me). I modified the old extrasites.py slightly to make the sequence uppercase. Rather than mess with that script any more, I write the results to disk as my.results.txt.

$ python extrasites.py > my.results.txt

So now the problem is to load to the two different textfiles with results and compare the data. A classic everyday bioinformatics problem. My output looks like this:

codon 123 GCG => GCT
AAAAAA GCG TTTGCG
AAAAAA GCT TTTGCG
HindIII AAGCTT

The original sequence is on the second line and the mutated sequence on the third. The affected codon is set off from the surrounding sequence by a space on each side.

The code to analyze the differences is pretty ugly. I listed it below but I hope you don't look at it unless you're stuck. The output, shown next, reveals that for the most part the two sets of results are congruent. EMBOSS results are output on a single line starting with 'E', while my results are output on three lines, the first starting with 'M'. The results were sorted by order of position in the sequence. The EMBOSS Base-Posn has been converted to a codon for consistency.

A few significant differences can be seen. Mostly they involve a single enzyme, KasI. I haven't sorted that out yet. Overall, I think the agreement is very good.

output:

$ python analyze.py 
E (40, 13, 'GGCGCC', 'KasI')
M ---  13  KasI
  AAAACG GCA CCTGTA
  AAAACG GCG CCTGTA
E (145, 48, 'GTCGAC', 'SalI')
M ---  48  SalI
  GTGCTG TCA ACCTGT
  GTGCTG TCG ACCTGT
E (199, 66, 'CTGCAG', 'PstI')
M ---  66  PstI
  AACCTG CAA GAAGCG
  AACCTG CAG GAAGCG
E (322, 107, 'TCTAGA', 'XbaI')
M ---  107  XbaI
  AGCGGT CTG GATTCA
  AGCGGT CTA GATTCA
E (331, 110, 'ACTAGT', 'SpeI')
M ---  110  SpeI
  GATTCA CTG GTGCTG
  GATTCA CTA GTGCTG
E (373, 124, 'AAGCTT', 'HindIII')
M ---  124  HindIII
  AAAAAA GCG TTTGCG
  AAAAAA GCT TTTGCG
E (481, 160, 'GGCGCC', 'KasI')
M ---  160  KasI
  ATCGGC GCT AGCGCC
  ATCGGC GCC AGCGCC
E (526, 175, 'AGATCT', 'BglII')
M ---  175  BglII
  GCCCGC CAA ATCTTT
  GCCCGC CAG ATCTTT
E (541, 180, 'GTCGAC', 'SalI')
M ---  180  SalI
  GAATCG CTC TCGACG
  GAATCG CTG TCGACG
M ---  180  SalI
  GAATCG CTC TCGACG
  GAATCG TTG TCGACG
E (571, 190, 'GGCGCC', 'KasI')
E (589, 196, 'ACTAGT', 'SpeI')
M ---  196  SpeI
  ATTGAA CTG GTGGCG
  ATTGAA CTA GTGGCG
E (685, 228, 'GGCGCC', 'KasI')
E (719, 239, 'CTGCAG', 'PstI')
M ---  240  PstI
  GCCCGT TTG CAGGAT
  GCCCGT CTG CAGGAT
M ---  248  SalI
  ATTATC AGT TCGACC
  ATTATC TCG TCGACC
M ---  295  MluI
  GCGAAC GCT TATCTT
  GCGAAC GCG TATCTT
E (904, 301, 'GTCGAC', 'SalI')
M ---  301  SalI
  AGCGTC GAT GATTTA
  AGCGTC GAC GATTTA
M ---  303  PstI
  GATGAT TTA CAGAGC
  GATGAT CTG CAGAGC
E (952, 317, 'CTGCAG', 'PstI')
E (952, 317, 'CTGCAG', 'PstI')
M ---  317  PstI
  CAGGCT GCG GCAGTA
  CAGGCT GCA GCAGTA
M ---  317  PstI
  CAGGCT GCG GCAGTA
  CAGGCT GCT GCAGTA
E (1021, 340, 'GGCGCC', 'KasI')
M ---  340  KasI
  GCCCAG GGG GCCAGC
  GCCCAG GGC GCCAGC
E (1130, 376, 'CTGCAG', 'PstI')
M ---  377  PstI
  GCCATC TTG CAGGAT
  GCCATC CTG CAGGAT
M ---  377  PstI
  GCCATC TTG CAGGAT
  GCCATC CTG CAGGAT
E (1135, 378, 'AGATCT')
M ---  378  BglII
  ATCTTG CAG GATCTG
  ATCTTG CAA GATCTG

code listing:

import REnzymes
RE = REnzymes.REnzymes()
good_enzymes = ['PstI', 'XbaI', 'SpeI','HindIII',
                'SalI', 'KasI', 'SpeI', 'PstI',
                 'BglII', 'MluI','Bcl']

def load_data(fn):
    FH = open(fn,'r')
    data = FH.read()
    FH.close()
    return data.strip()
#===================================
# part 1:  EMBOSS results
data = load_data('emboss.results.txt')
data = data.split('\n\n')[2]
assert data[:7] == '  Start'
eL = list()
for e in data.strip().split('\n')[1:]:
    e = e.split()
    x, j, flag, enz, seq, i = e[:6]
    if not enz in good_enzymes:  continue
    if not '+' in flag:  continue
    codon = int(i)/3
    i = int(i)+1
    eL.append((i,codon, seq, enz))
#===================================
# part 2:  my results
data = load_data('my.results.txt')
data = data.split('\n\n')[:-1]
mL = list()
for entry in data:
    lines = entry.split('\n')
    codon = lines[0].split()[1]
    # bug in original
    codon = int(codon) + 1
    enz,seq = lines[3].split()
    if not enz in good_enzymes:
        continue
    L = ['---', str(codon), enz]
    L += [lines[1].strip()]
    L += [lines[2].strip()]
    mL.append(L)
#===================================
# part 3:  show
e = eL.pop(0)
m = mL.pop(0)
while eL or mL:
    if int(m[1]) < int(e[1]):
        print 'M', '  '.join(m[:3])
        print ' ', '\n  '.join(m[3:])
        m = mL.pop(0)
    else:
        print 'E', e
        e = eL.pop(0)
    if not mL:
        print 'E', e[:3]
    if not eL:
        print 'M', '  '.join(m[:3])
        print ' ', '\n  '.join(m[3:])

Friday, November 12, 2010

Code to find restriction sites

Having downloaded a list of restriction enzymes from REBASE earlier today for another post, I couldn't resist writing a short program to find restriction sites in a DNA sequence. The input enzyme data looks like this:

AarI                           CACCTGC (4/8)
AatII                          GACGT^C
AbsI                           CC^TCGAGG
AccI                           GT^MKAC
AceIII                         CAGCTC (7/11)
AciI                           CCGC (-3/-1)

We ignore the numbers (4/8) etc. To make it more interesting, I used a simple form of Python's regular expressions (docs here), which are pre-compiled.

The results can be limited to non-degenerate six-cutters (or longer) and sorted by the index of the match. It took a bit more than an hour, and was great fun. Naturally, the code hasn't been tested carefully. (Note: 0-based indexing).

output:

$ python script.py
Eco47III   AGCGCT             11  AGCGCT
HaeII      RGCGCY             11  AGCGCT
TsoI       TARCCA             23  TAACCA
HgiCI      GGYRCC             35  GGCACC
AflIII     ACRYGT             56  ACGCGT
MluI       ACGCGT             56  ACGCGT
AclI       AACGTT             62  AACGTT
BclI       TGATCA             83  TGATCA
BspGI      CTGGAC             93  CTGGAC
MstI       TGCGCA            107  TGCGCA
BsgI       GTGCAG            120  GTGCAG
HindII     GTYRAC            140  GTCAAC
BspMI      ACCTGC            190  ACCTGC

code listing:

import re
sample_dna = '''
 ATGACCCTTTTAGCGCTCGGTATTAACCATAAAACGGCACCTGTATCGCT
 GCGAGAACGCGTAACGTTTTCGCCGGACACGCTTGATCAGGCGCTGGACA
 GCCTGCTTGCGCAGCCAATGGTGCAGGGCGGGGTCGTGCTGTCAACCTGT
 AACCGTACAGAGCTGTATCTGAGCGTGGAAGAGCAGGATAACCTGCAAGA'''

# load data from downloaded file
# http://rebase.neb.com/rebase/link_proto
def load_data(fn = 'link_proto.txt'):
    FH = open(fn,'r')
    data = FH.read()
    FH.close()
    L = data.strip().split('\n\n')
    return L

# parse data into names and sites
def preprocess(data):
    L = data.strip().split('\n')
    names = list()
    sites = list()
    for e in L:
        n,s = e.split(' ',1)
        names.append(n)
        words = s.strip().split()
        if words[0][0] == '(':
            sites.append(words[1])
        else:
            sites.append(words[0])
    return names,sites

# codes for degenerate positions
# http://www.bioinformatics.org/sms/iupac.html 
def get_pattern(s):
    D = { 'A':'A','C':'C','G':'G','T':'T',
          'R':'[AG]','Y':'[CT]','N':'.',
          'S':'[GC]','W':'[AT]','K':'[GT]',
          'M':'[AC]','B':'[CGT]','D':'[AGT]',
          'H':'[ACT]','V':'[ACG]',
          '^':''}
    rL = [D[c] for c in s]
    return ''.join(rL)

# dictionary of pre-compiled regexps
def make_dict(data):
    names,sites = preprocess(data)
    D = dict()
    for n,s in zip(names,sites):
        p = get_pattern(s)
        p = re.compile(p)
        i = s.find('^')
        s = s.replace('^','')
        rD = { 'name':n,'site':s,
               'pattern':p,'i':i }
        D[n] = rD
    return D
    
def search(dna,D,minlength=4,ignore_ambig=True):
    N = max([len(n) for n in D.keys()])
    M = max([len(D[n]['site']) for n in D])
    rL = list()
    for n in D:
        rD = D[n]
        n,p,s = rD['name'], rD['pattern'],rD['site']
        if len(s) < minlength:  
            continue
        if ignore_ambig and 'N' in s:
            continue
        m = p.search(dna)
        if m:
            i = m.start()
            j = i + len(rD['site'])
            e = [n.ljust(N),s.ljust(M)]
            e += [i,dna[i:j]]
            rL.append(e)
        rL.sort()
    return rL

def show(result):
    def f(s):  return s[2]   # index of match
    result = sorted(result, key=f)
    for line in result:
        # left i as int, so convert to str
        i = line[2]
        line[2] = str(i).rjust(4)
        print '  '.join(line)
            
if __name__ == '__main__':
    dna = ''.join([c for c in sample_dna if c in 'ACGT'])
    data = load_data()
    L = data[3].strip()  # type II enzymes only
    D = make_dict(L)
    
    result = search(dna,D,minlength=6)
    show(result)