Installing libsvm-3.0 for Python on OSX 10.6

This was very frustrating to solve and I eventually had to get my friend Kieran (thanks Kieran!) to help me.  Basically the vanilla installation instructions that come with are insufficient (at least for me) to get a working module.  Here are the steps that were required to get everything to work.

1. Run 'make' in the libsvm-3.0 directory

2. Run 'make' in the libsvm-3.0/python directory

In the libsvm-3.0 directory there now should be a .so2 file

3. Create a new directory in your site-packages directory (your pythonpath) called libsvm

4. Copy the .so2 file from libsvm-3.0 and the svm.py, svm.pyc, svmutil.py files from libsvm-3.0/python to site-packages/libsvm

There are a couple things that are missing so now we need to make them.

5. In site-packages/libsvm create a file called __init__.py.  This is an empty file, but it is necessary to get the directory recognized as a python module.

6. Edit svm.py and add the following two lines after the other import statements at the top of the file:

          import os.path
    _PATH = os.path.join( *os.path.split(__file__)[:-1] )


7.  At around line 7 you will see this statement


          # For unix the prefix 'lib' is not considered.

    if find_library('svm'):
       libsvm = CDLL(find_library('svm'))
    elif find_library('libsvm'):
       libsvm = CDLL(find_library('libsvm'))
    else:
       if sys.platform == 'win32':
          libsvm = CDLL('../windows/libsvm.dll')
       else:
          libsvm = CDLL('../libsvm.so.2')

8.  Change this to look like this:


    # For unix the prefix 'lib' is not considered.
    if find_library('svm'):
  libsvm = CDLL(find_library('svm'))
    elif find_library('libsvm'):
  libsvm = CDLL(find_library('libsvm'))
    else:
  if sys.platform == 'win32':
  libsvm = CDLL(os.path.join(_PATH,'windows','libsvm.dll'))
  else:
  libsvm = CDLL(os.path.join(_PATH,'libsvm.so.2'))

9. Once you save svm.py, you should be able to fire up a python interpreter and do 'from libsvm import svm'.  If that works, and a dir(svm) shows you a ton of functions, then you are good to go.

Generate all possible proteins from ambiguous DNA

This had me stumped for awhile, but this works pretty well.  Does NOT handle stop codons or gap characters like '-'.  Requires BioPython

import itertools
from Bio.Seq import Seq
from Bio.Data import CodonTable
from Bio.Data import IUPACData</pre>

# Takes Bio.Seq.Seq object as input
# Returns list of all possible proteins
# Assumes sequence is in frame +1
def generateProtFromAmbiguousDNA(s):
   std_nt = CodonTable.unambiguous_dna_by_name["Standard"]
   nonstd = IUPACData.ambiguous_dna_values
   aa_trans = []
   for i in range(0,len(s),3):
      codon = s.tostring()[i:i+3]
      aa = CodonTable.list_possible_proteins(codon,std_nt.forward_table,nonstd) 
      aa_trans.append(aa)
   proteins = list(itertools.product(*aa_trans))
   possible_proteins = []
   for x in proteins:
      possible_proteins.append("".join(x))
   return possible_proteins

def main():
   a = Seq('ATGGCARTTGTAHAC')
   print "DNA: ",a.tostring()
   print "Proteins:"
   foo = generateProtFromAmbiguousDNA(a)
   for s in foo: print s

if __name__ == '__main__':
   main()

Creating a quick codon table

I didn't think this up, the code comes from Peter Collingridge here.  But it is rather elegant.

bases = ['t', 'c', 'a', 'g']
codons = [a+b+c for a in bases for b in bases for c in bases]
amino_acids = "F F L L S S S S Y Y stop stop C C stop W L L L L P P P P H H Q Q R R R R I I I M T T T T N N K K S S R R V V V V A A A A D D E E G G G G".split(' ')
codon_table = dict(zip(codons, amino_acids))

Update the locate database on the Mac

This is the command for updating the locate database on the OSX system.

sudo /usr/libexec/locate.updatedb


I should figure out how to make this run everyday.

Connecting to PostgreSQL with Python and Psycopg2

Basic syntax for making a database connection, executing and retrieving data:

import psycopg2 as pg

# create database connection
try:
   conn = pg.connect("dbname='template1' user='dbuser' host='localhost' password='dbpass'")
except:
   print "Unable to connect to database"


# create database cursor
cur = conn.cursor()


# execute SQL and fetch results
cur.execute("""SELECT datname from pg_database""")
rows = cur.fetchall()


print "\nShow database results:\n"
for row in rows:
   print row[0]

Density plots in R

To plot the distribution of scores in R using the density function, use

plot(density(dat$scores))

To overlay scores for another variable onto the first plot, use

lines(density(dat$other_scores), col="red")

Simple enough!

Using Multiprocessing in Python

For a newbie to multi CPU processing, I have to say that the python 2.7 documentation is difficult to understand at best and incomplete at worst.  Searching the web I came across two tutorials that I found to be very helpful.

The first is by Doug Hellmann in his Python Module of the Week series (PyMOTW), he gave much clearer examples of how the multiprocessing module worked.

However, I still needed more information.  I found enough to get me over the hump in solving my problem from Norman Matloff's (pdf link) tutorial from UC Davis, called "Programming on Parallel Machines".  Chapter 3 is called "The Python Threads and Multiprocessing Modules".  Not exhaustive, but very helpful.  I will probably be referring to it again as I move the code from a single multiprocessor machine to a cluster.

Other information on other types of Parallel Processing (Cluster, Cloud, Grid) and the related python libraries can be found here at the PythonWiki