Plotting data with matplotlib

In continuing my slow migration away from “office-like” tools for working with my data, I’ve been taking a look lately at matplotlib. I’ve banged together a rough script to do some simple data plotting with a bit of flexibility:

#! /usr/bin/env python
# http://biostumblematic.wordpress.com

# An interface to matplotlib

# Import modules
import csv, sys
import matplotlib.pyplot as plt
import numpy as np

# Introduce the program
print '-'*60
print 'Your data should be in CSV format, with Y-values'
print 'in odd columns and X-values in even columns.'
print 'If your file contains a header row, these will be'
print 'automatically detected'
print '-'*60

# Open the data
datafile = sys.argv[1]
f = open(datafile, 'r')

# Check to see if the file starts with headers or data:
dialect = csv.Sniffer().has_header(f.read(1024))
f.seek(0)
reader = csv.reader(f)

# Assign the data to series via a dict
if dialect is True:
	reader.next() # Move down a line to skip headers
else:
	pass

series_dict = {}
for row in reader:
	i = 0
	for column in row:
		i += 1
		if series_dict.has_key(i):
			try:
				series_dict[i].append(float(column))
			except ValueError:
				pass
		else:
			series_dict[i] = [float(column)]
# Plot each data series
num_cols = len(series_dict)
i = 1
while i < num_cols:
	plt.plot(series_dict[i], series_dict[i+1], 'o')
	i += 2 

# Get axis labels
xaxis_label = raw_input('X-axis label > ')
yaxis_label = raw_input('Y-axis label > ')

# Show the plot
plt.ylabel(yaxis_label)
plt.xlabel(xaxis_label)
plt.show()

# Enter loop for customizing appearance

# Stop
f.close()

As-is this will read in a CSV file of any number of columns and plot them as Y values/X values (alternating).

Some things that feel nasty:

• Having to use the dictionaries to get the column data assembled. I feel like the CSV reader module should have a “transpose” function
• The section near the end where I’m generating the different plots by iterating over the number of columns.

Some things that would be nice to implement, but I haven’t figured out yet:

• More differentiation of the appearance for each series’ plot
• Automatic generation of a legend using headers for the X-values from the initial file (or else requested from the user at run-time if not in the file)

Using conky + python to display twitter feeds

No posts for a while because I haven’t actually been writing anything new. Biopython has solved many of my day-to-day problems, and I’m in love with SeqIO.

Today is Canada Day and it’s pretty quiet around the lab, so I thought I’d try to write something that would let me do two things:
First, I want to be able to view my twitter feed using conky, and secondly I’d like to be able to send updates from the console.

This also gives me the chance to work on some fundamentals like interfacing with APIs and passing options to scripts from the command line. There, I totally justified it!

To start off, I installed python-setuptools (from the Ubuntu repo), simplejson and the python-twitter API interface. To install the last two you just download the archives, extract them, and then run the following two commands from within their folders:

python setup.py build
sudo python setup.py install

Let’s start off with a pretty basic framework. This script should print out your latest five friends’ updates to the terminal, and the struts are in place to eventually add post capability:

#! /usr/bin/env python
# http://biostumblematic.wordpress.com
# Simple twitter interface

# Change the following two lines with your credentials
user = 'username'
pw = 'password'

num_statuses = 5 # Changes number of statuses to show

import sys, twitter
api = twitter.Api(username=user, password=pw)

if sys.argv[1] == '-l':
    timeline = api.GetFriendsTimeline(user)
    i=0
    while i < num_statuses:
        print timeline[i].user.name
        print timeline[i].text
        print '\n'
        i+=1

elif sys.argv[1] == '-p':
    pass

else:
    print 'Invalid input'
    print 'Allowed options are:'
    print '-p (to post an update)'
    print '-l (to list friend statuses)'
    sys.exit(2)

Adding the update functionality is facile. Just change the code as follows:

elif sys.argv[1] == '-p':
    status = api.PostUpdate(sys.argv[2])
    print 'Twitter status updated'

The only caveat in doing it this way is that the status update entered at the command line must be passed as a string, with quotation marks around it. Otherwise this will post a one word update, which is terse even by twitter standards.

This is already fully functional in the terminal, so the last step is dumping out the statuses to a file which conky can read. Here’s the code I used to make the text file:

if sys.argv[1] == '-l':
    timeline = api.GetFriendsTimeline(user)
    i=0
    output = open(os.environ['HOME']+'/tweets.txt', 'w')
    while i < num_statuses:
        output.write(timeline[i].user.name+'\n')
        output.write(timeline[i].text+'\n')
        output.write('\n')
        i+=1
    output.close()

We then have conky read it using a file like this (I named in .conkytweets and placed it in my home directory, make sure to change your home directory below):

use_xft yes
xftfont MyriadPro-Regular:size=8
alignment top_left
xftalpha 0.8
own_window yes
own_window_type override
own_window_transparent yes
own_window_hints undecorated,below,sticky,skip_taskbar,skip_pager
double_buffer yes
draw_shades no
draw_outline no
draw_borders no
stippled_borders 10
border_margin 4
border_width 0
default_shade_color black
default_outline_color black
use_spacer right
no_buffers no
uppercase no
default_color 222222
maximum_width 200
minimum_size 200 5
gap_y 400
gap_x 10
text_buffer_size 1024

TEXT
${font size=9}Latest Tweets:
${color}${font}${execi 600 cat /home/jason/tweets.txt | fold -w 35}

Here is a screenshot of the output on my monitor (sorry for the blur over the tasks, there are some research details in there I just didn’t feel like posting for the whole world atm)

demo of conkytweets script

There are a few things that don’t work very well. To my knowledge, you can’t include clickable links in conky, so URLs in tweets don’t do anything. The textwrap in conky is also a bit wonky, but I don’t know that there is a nice fix for that. I suppose one option would be to modify the text file that the twitter script generates, but I’ll leave that as an exercise to the reader.

The simplest way to use this is to add a link to your path. For me it was:

cd /usr/bin/
sudo ln -s ~/scripts/pytwit.py pytwit

Then you can use it from anywhere with either:

pytwit -p 'My awesome twitter post'

or

pytwit -l

For extra points, you can add the listing to your crontab as follows:

sudo gedit /etc/crontab
*/15 *	* * *	jason pytwit -l

This will update the statuses every 15 minutes.

The full script is available on github, and I welcome any additions/modifications/improvements as always.

Measuring identities of aligned protein sequences with BioPython

For some reason it seems that every program which will output a percentage of the identity between two proteins will also align them itself – therefore screwing up any alignment which you’ve already made. I knocked up a short script over the weekend which will read in a FASTA-formatted alignment and output the percent identity of all of the proteins in it to the first one in the file.

I couldn’t find a built-in way to do this all in BioPython, but I did use it to parse the seqences out of the alignment. The rest of the work is just brute force string crunching.

#!/usr/bin/env python
# http://biostumblematic.wordpress.com

import string
from Bio import AlignIO

# change input.fasta to match your alignment
input_handle = open("input.fasta", "rU")
alignment = AlignIO.read(input_handle, "fasta")

j=0 # counts positions in first sequence
i=0 # counts identity hits
for record in alignment:
    for amino_acid in record.seq:
        if amino_acid == '-':
            pass
        else:
            if amino_acid == alignment[0].seq[j]:
                i += 1
        j += 1
    j = 0
    seq = str(record.seq)
    gap_strip = seq.replace('-', '')
    percent = 100*i/len(gap_strip)
    print record.id+' '+str(percent)
    i=0

I didn’t implement similarity here, but it gets the basic job done. This script is available on GitHub as seqhomology.py

A worked example in BioPython: From cDNA to protein and back again

I’m not sure if this behavior is normal, but I find that I learn a new system best by choosing a real-life problem that I need to solve and applying the new method in order to solve it. This inevitably means that I’ll probably be doing things in a non-efficient way (since I’m a noobie), but code can always be refined later.

Here is the problem I have in front of me today: I have a series of proteins from which I’d like to isolate (via cloning) a certain domain. The cDNA clones of the full length proteins are available from the IMAGE consortium. Unfortunately these aren’t completely “clean” cDNAs; there tends to be some extraneous sequence on both ends of the gene.

The plan of action goes something like this:
The starting materials are the cDNA sequence, the amino acid sequence of the protein, and the residue ranges of the domain of interest. So what I’d like to do is to check each frame of the cDNA to find the one matching the translated protein sequence, then extract just the cDNA coding for the domain I’d like to isolate. I can then (independently) design PCR primers for this domain.

You’re probably thinking that this could be done manually (and of course that’s true), but I find this painstaking work. Also it gives me a chance to play around with the SeqIO functions of BioPython a bit

Enough introduction, let’s get to work. The protein I’ll be using for this exercise is IKBKB. This is a 756 amino acid protein; I’ll be trying to get the cDNA for the protein kinase domain from residues 15-300. The IMAGE clone ID is 5784717.

Baby step 1 – find the ORF we’re interested in

#! /usr/bin/env python

# http://biostumblematic.wordpress.com

# Extraction of the cDNA
# for a given protein domain

import re
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC

input_cdna = raw_input('Paste your cDNA sequence >> ')
input_search = raw_input('What are the first amino acids of the protein? >> ')

cdna = Seq(input_cdna, IUPAC.unambiguous_dna)

i=0
while i < 3:
    frame = cdna[i:150]
    trans = frame.translate()
    orf_find = re.search(input_search, str(trans))
    if orf_find:
        orf_frame = i+1
    else:
        pass
    i += 1
print 'The protein is coded in frame '+str(orf_frame)

Given the input of the cDNA and the first 4 residues (MSWS), this outputs the right answer:

The protein is coded in frame 3

Note that I’m only checking the first 50 residues of the cDNA (see line 19) Hopefully this is enough to catch the protein of interest (it would be a lot of extraneous 5′ sequence if not).

Obviously this is not enough sexy for Biopython. It’s cumbersome to have to type in the starting sequence of the protein you’re interested in, so why don’t we let SeqIO handle that for us via the SwissProt code? Also, we’ll change a couple of things to enable automation of the full list later.

First, I made a file called ‘test.csv’ which has a single line consisting of SwissProt ID,cDNA:

O14920,atagccccggg[...]

Then I modified the script like so:

import re, csv
from Bio import ExPASy, SeqIO
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC

reader = csv.reader(open('test.csv'))
for row in reader:
    input_prot = row[0]
    get_prot = ExPASy.get_sprot_raw(input_prot)
    prot_obj = SeqIO.read(get_prot, "swiss")
    get_prot.close()
    prot_seq = prot_obj.seq
    prot_start = prot_seq[0:4]

    cdna = Seq(row[1], IUPAC.unambiguous_dna)

    i=0
    while i < 3:
        frame = cdna[i:150]
        trans = frame.translate()
        orf_find = re.search(str(prot_start), str(trans))
        if orf_find:
            orf_frame = i+1
        else:
            pass
        i += 1
    print 'The protein is coded in frame '+str(orf_frame)

Biopython grabs the protein sequence from the web using the SwissProt ID. The prot_start variable takes just the first few residues and uses that as the search term for the regular expression later on. Now there is no command line input, as everything is done via the CSV file. This will iterate over lines in the CSV file to do multiple proteins. Right now, however, we would just get a long list of “The protein is coded in frame X” lines, which is less than useful. Time to take care of that.

In this case the domain is annotated in SwissProt already. This means that I could use the built-in parsing function of BioPython to select the domain, however I have some custom annotations for other proteins in my list that make this not a good idea in this case. Instead let’s just make some minor modifications to our input CSV and existing script. The new CSV includes the start and stop residues of interest:

O14920,15,300,atagccccgggttt[...]

Now I just modify the top of the script to take into account the new structure of the CSV:

reader = csv.reader(open('test.csv'))
for row in reader:
    input_prot = row[0]
    get_prot = ExPASy.get_sprot_raw(input_prot)
    prot_obj = SeqIO.read(get_prot, "swiss")
    get_prot.close()
    prot_seq = prot_obj.seq
    prot_domain = prot_seq[int(row[1])-1:int(row[2])]
    cdna = Seq(row[3], IUPAC.unambiguous_dna)

and adjust what happens if the script finds a match to the domain sequence:

        if orf_find:
            trans_split = re.split('('+str(prot_domain)+')', str(trans))
            cdna_start = len(trans_split[0])*3
            cdna_stop = cdna_start + len(trans_split[1])*3
            cdna_extracted = frame[cdna_start:cdna_stop]
            print cdna_extracted

And this gets the job done! This prints a cDNA sequence which, when translated back, matches the domain of interest.

The last part there feels sloppy. All I’m doing is counting the number of amino acids that come out of the translation before the start of the domain, then multiplying by three, and getting the cDNA start from this. I feel like there must should be a way to transition more effectively between protein and cDNA sequence.

The entire script, slightly modified to write out the results to a new CSV file, is available over on GitHub. I hope you found the post interesting and look forward to your comments.

Why the #^@% have I not done more with Biopython before now?

Yeah, I’ve heard of it. Biopython is A python module package (thanks Chris) that’s written to help with doing computational biology. To my utter dismay I somewhat ignored it, being the “ll just brew it myself” type. What a mistake.

Today I was trying to wrangle some DNA and protein sequences and realized that this might be something covered by Biopython. It’s even better than that. You want tasty yum yums? How about a reverse complementer in 3 lines of code? I even formatted it so it looks nice on the terminal:

from Bio.Seq import Seq
sequence = Seq(raw_input('Paste your DNA sequence >> '))
print '\nReverse Complement\n------------------\n'+sequence.reverse_complement()

The very next bit of code in the tutorial replaced a ~50 line program I had cobbled together (and which still wasn’t working exactly the way I wanted) into this beauty:

#! /usr/bin/env python

# Biopython can automatically parse FASTA
# as well as many other "standard" biological formats

from Bio import SeqIO
inputfile = open('myproteins_fasta.txt')

for seq_record in SeqIO.parse(inputfile, 'fasta'):
    print seq_record.id
    print repr(seq_record.seq)
    print len(seq_record)
inputfile.close

BOOM, FASTA reader.

I’m just getting started on reading the documentation, but so far I’m really impressed (and not a little bit sheepish at my previous obstinance). Expect to see some Biopython examples in the coming days

Parsing UniProt XML with ElementTree

I’ve written about ElementTree before, and it really is a handy tool. I took the output NCBI GI numbers from my previous post and used them in concert with the ID mapper at UniProt to get a listing of the proteins. UniProt kindly allows you to download this subset in XML, which I did in order to quickly extract the information I was interested in.

Here’s what a bit of the XML looks like:

Some throwaway code for playing with poorly formatted text files

I managed to find a paper in which some of the analysis I’ve been working on had been done. Unfortunately the raw results of the analysis were just that – raw. Specifically they had been dumped into a 6.8 MB text file as a supplement to the paper.

In order to extract the information I was interested in, and to prove to people who read this that I don’t solve all of my problems with Python, I thought I’d share the quick code I used.

First of all, I wanted all of the lines that reported proteins from humans. This turned out to be workable by running:

cat infile.txt | grep 'Homo sapiens' > oufile.txt

This gave me a long list which helpfully had each line starting with the NCBI GI number for the protein of interest. To extract the GI numbers alone involved:

cat outfile.txt | cut -c 1-11 > GI_list.txt

then to trim the whitespace:

sed 's/^[ \t]*//;s/[ \t]*$//' GI_list.txt > GI_list.txt

(this last one took some help from the handy sed one-liners page)

The entire process took about 1/4 of the time I’ve just taken writing it up, and I now have a nicely-formatted 11 kb file which I can use as input to my next round of tasks.

“Inheritance” partially solved

I still don’t claim to understand inheritance well (although it should be simple), but I think that part of my trouble was stemming from differences between Django and Python. Regardless, I managed to solve my problem from the last post with something like this:

class Protein(models.Model):
    ...
class Domain(models.Model):
    protein = models.ForeignKey(Protein)
    ...
    def domain_percent(self):
        m = re.match('^RING', self.domain)
        if m:
            p = Protein.objects.get(pk=self.protein_id)
            p_clean = re.sub("\s+", "", p.sequence)
            p_len = float(len(p_clean))
            domain_l = float(self.domain_stop - self.domain_start)
            perc = int(100*(domain_l/p_len))
            return perc
        else:
            pass

the key line is

p = Protein.objects.get(pk=self.protein_id)

. This says “get all of the objects from the related class “Protein” which have a primary key (pk) equal to the protein_id of this domain. Whew.

As usual, the key information I needed was in the part of the tutorial that I skimmed. I’ve spent most of the day today struggling with querysets, but I’m still a long way from maximizing them as well. I’ve managed to get things working “well enough” for the time being – I need to stop tinkering with the backend of this for now and actually populate/extract data if I want to be able to wow my advisor at our individual meeting tomorrow.

I fundamentally don’t grok inheritance in python

<insert ramble about how long it’s been since the last post here>

So I’m back to working with Django, and once again something that seems simple is killing me.  I simply want to be able to define the various conserved domains in the protein, then be able to automatically pull out the amino acid sequence of a given domain.  The problem (which I’ve dealt with before) is getting the python classes to talk to one another.

Here is some simple code to illustrate:

class Protein(models.Model):
    uniprot = models.CharField('Uniprot accession code', max_length=6, blank=True)
    sequence = models.TextField('Amino acid sequence')

class Domain(models.Model):
    protein = models.ForeignKey(Protein)
    domain = models.CharField('Domain Name', max_length=25, blank=True)
    domain_start = models.IntegerField('Domain start residue', blank=True)
    domain_stop = models.IntegerField('Domain end residue', blank=True)

So what I want to do looks like this in my head, but doesn’t work:

class Domain(models.Model):
...
    def domain_FASTA(self):
        amino_clean = re.sub("\s+", "", protein.sequence)
        domain_sstart = self.domain_start - 1 # lists start from 0
        domain_sstop = self.domain_stop - 1
        domain_seq = amino_clean[domain_sstart:domain_sstop]
        fasta = '> '+protein.uniprot+' - '+self.domain+' '+str(self.domain_start)+'-'+str(self.domain_stop)+'\n'+domain_seq
        return fasta

I’ve been banging on this for hours, and have managed to fairly thoroughly break the things that I did have working. Fortunately I’m older and wiser than I once was, so I managed to restore them. I’m pretty well stumped on how to solve this one, though.

4 Hour Liveblog: Website development with Django and Pinax

One of the projects I’ll be involved in at the start of my new post-doc is developing a website. This site will be a portal to enable researchers in our field to communicate, as well as a gateway to the database I will (likely) be developing.

I’ve tinkered with Django before and really liked it. I’ve recently heard about Pinax, which is an open-source set of Django applications that are designed for general “plug and play”-like use on a variety of websites. Many of the features of Pinax are elements that would be great to implement on this portal site I need to work on.

I’ve heard quite a bit about how fast it is to develop sites using these methods, so I’ve decided to take 4 hours out of my Saturday and put this to the test. Of course I’m a super-green programmer, so it’s a tough test indeed. I’ll be updating this post as I go with progress reports, so stay tuned!

10:30-10:45 AM – Setting up the environment
I’m working on two machines today. My linux box is where the actual Django/Pinax will live, and I’m using a Windows machine to update the blog and do graphics work (I’m still a slave to Photoshop).

I’ve already got a recent SVN version of Django installed, so the first real step is installing Pinax. First let’s download the SVN version:

svn checkout http://svn.pinaxproject.com/pinax/trunk/ pinax

The installation docs for Pinax include some simple instructions for obtaining dependencies. Following the instructions verbatim gives you both a recent copy of Django as well as Pinax itself, so you can probable just start from there if you are new to the system.

You can then use a couple of methods to install the various dependencies. Since I’m doing this the quick way, I just ran:

pip install --requirement src/pinax/requirements/external_apps.txt

Which should install everything. This cranked away in the terminal for a bit.

After it finished, I changed into the src/pinax/pinax/projects/complete_project/ directory and ran:

./manage.py syncdb

and

./manage.py runserver

and here is what we see:

Vanilla pinax installation

Not bad at all for 15 minutes of work. Next up: investigating the Pinax back end and customization.

10:45 – 11:15: Exploration
During the installation process, I was asked to create a super user. These credentials can be used to log in to the vanilla pinax site. On the login screen, we already see one of the features of pinax – integration with OpenID:

pinax login screen

Once logged in, you are directed to a page which gives you some user functions. What I want to see first however is the admin area, which is available as a link at the top right of the screen. This is just the baked-in admin that Django itself includes (which is really great, don’t get me wrong). You can click on the various modules to change them, add new items, etc. To be honest, this isn’t as exciting as I thought it might be; let’s go back to the user page.

Just exploring the tabs along the top, it’s very simple to edit your user profile, create projects, etc.

Twitter clone in Pinax

I did run into a few bugs here. When trying to add a bookmark, I got an error which broke the page, giving a template syntax error:

Caught an exception while rendering: 'NoneType' object has no attribute 'append'

The “Locations” tab also coughed up an error, but this is simply because I hadn’t entered an API key for Yahoo maps yet.

When trying to add events to a tribe’s calendar, it initially didn’t work at all. I had to first enter an event in the admin screen to get the calendar link on the user side to work at all. Once this was in place, trying to add an event was difficult because the date/time fields were just unguided text entry that would give an error stating “Enter a valid date/time.” whenever I would try to put in the information.

Looks like we will have to do some troubleshooting.

11:15 – 11:45: Moving out of the nest
Before we get into this, however, let’s get ourselves into a better position from which to work from.
I copied the pinax “basic project” folder out into a new directory. I edited the settings.py file to point at the Pinax installation, but got an import error when trying to view this new project:

ImportError at /

No module named basic_project.urls

Next I tried going into my pinax bin directory and using the following command:

python pinax-clone-project.py ../projects/basic_project /home/user/path/to/myproject/

I then changed into the new directory and edited the settings.py file (line 83) to read:

ROOT_URLCONF = 'myproject.urls'

This worked, and I could now log in to the basic site. Next I’ll re-enable some of the apps.

Hour 1 is done

11:45 – 12:30: Enabling applications and customizing the template
Since I cloned the basic project, I want to turn back on some of the applications. I opened up the settings.py files for the “complete project” included with Pinax and that for my project, and pasted in the applications I wanted under INSTALLED_APPS. I also had to copy some code over into TEMPLATE_LOADERS, MIDDLEWARE_CLASSES, TEMPLATE_CONTEXT_PROCESSORS and COMBINED_INBOX_COUNT_SOURCES. In the end the changes were so extensive that I just gave in and copied over the entire text of the complete_project settings.py file and commented out the stuff I didn’t want.

Unfortunately, when I would try to run syncdb, I’d get an error:

Error: No module named photos

which I couldn’t trace back by commenting out certain lines, etc. At this point I’m starting to get a hair frustrated.

Back to square one. I deleted that project and just cloned over the complete project as described above, and we were back in business.

Now let’s try to get rid of the applications I don’t want: swaps and locations. I had already commented these out under INSTALLED_APPS in settings.py, but they still show up in the template. This should be an easy fix… let’s just open up site_base.html in the templates directory and delete the code for these tabs.

We don’t even have to restart the development server for this one; just refresh the page and they are gone!

That took much longer than I had hoped… I need a short break.

12:30 – 2:30: Theming and visual customization
Rather than break this out into small chunks, I’m just going to lump the second half under theming and customization. Obviously I don’t want the vanilla look complete with Pinax logo at the top, so let’s dive into the Django templating system.

Firstly, I made a directory in the myproject folder called site_media, and placed a custom logo in a subfolder of this called images. To display this image, I just had to modify site_base.html under the templates directory to read:

{% block logo_link_image %}<a href="{% url home %}"><img src="{{ MEDIA_URL }}/images/logo.png" alt="" /></a>{% endblock %}

I also had to change a little bit of myproject/media/base.css in order to avoid resizing the new logo.

Most of the rest of my time was spent playing around with the base.css file, just getting a feel for what defined the different bits of the site. I think I got a little burnt out on some of the minor bugs, and just didn’t have the drive to get my hands dirty mucking about at the lower levels any more today.

That isn’t to say that I’m not very impressed with Pinax. I haven’t even mentioned some of the nice features (like built-in RSS feeds for the twitter clone). There are some things that I’d like to see added, like auto-population of slug fields, however I plan on working some more with this system over the coming weeks to see if I can’t iron out some of the kinks & problems. I’m just a bit of a slow learner, so it’s going to take me longer than this afternoon to sort out exactly what changes need to be made and where to make them.