Noel O'Blog

Chemistry in R

2008-07-23T08:47:00.000-07:00

For many cheminformaticians, R is the preferred way of analysing multivariate data and developing predictive models. However, it is not so widely known that there are R packages available that are directly aimed at handling chemical data.

Over the last few years, Rajarshi Guha (Indiana University) has been doing some nice work integrating the CDK and R. His publication in J. Stat. Soft., Chemical Informatics Functionality in R, describes the rcdk and rpubchem packages. The rcdk package allows the user to read SDF files directly into R, calculate fingerprints and descriptors, calculate Tanimoto values, view molecules in 2D (JChemPaint) and 3D (Jmol), calculate SMILES strings, and access the property fields of file formats. The rpubchem package is focused on downloading compounds, property values and assay data from PubChem. See also articles in R news and CDK news [1], [2], [3].

A more recent development is ChemmineR, described in the latest issue of Bioinformatics: "ChemmineR: a compound mining framework for R". The authors appear to be unaware of the earlier work by Rajarshi, and so there is no comparison of available features. However, based on the documentation on their website, it seems that much of the functionality revolves around a type of fingerprint called atom-pair descriptors (APD). SDF files, when read in, are converted to a database of APDs and these can be used for similarity searching, clustering, removal of duplicates and so on. Sets of molecules can be visualised using a web connection to the ChemMine portal (I'm not sure what software is used). According to the documentation, future work will include descriptor calculation with JOELib.

So, there you have it. An exhaustive survey of the two available methods for bringing chemistry into R. Is the time ripe for a cheminformatics equivalent to Bioconductor?

Review of "IronPython in Action"

2008-07-15T14:00:00.000-07:00

The three most well-known implementations of Python are the original C implementation (CPython), one in Java (Jython), and one in .NET (IronPython). Jython makes it easy to write Python programs that use Java classes and to test a Java library at an interactive prompt. Similarly, IronPython allows Python programs to be written that use the .NET libraries (these are the same libraries used by C# and Visual Basic).

IronPython (sometimes abbreviated as FePy) is relatively new and had its first release in 2006. It is Open Source and developed by Microsoft (by the same developer who started Jython). Although .NET is strictly for Windows, IronPython can also run on top of Mono, the open source implementation of .NET which is available cross-platform.

As a cheminformatician should you be interested in IronPython or .NET? Well, Dimitris Agrafiotis thinks so. And Eli Lilly recently open sourced a life science platform implemented in .NET (here's the SF website, but where's the mailing list?). I'm not yet convinced, but I'm keeping an open mind. I'm helping Matt Sprague to prepare C# bindings for OpenBabel. I think this will make OpenBabel the first cheminformatics library accessible from C# (although probably someone will contradict me below).

At this point, if you're still interested, you should check out the first chapter of "IronPython in Action", which you can read on the web. This explains in more detail the background to IronPython and why you might want to use it. IronPython in Action will be published later this year by Manning Press, and is written by Michael Foord and Christian Muirhead. Foord (aka Fuzzyman) is a very active blogger on Python and FePy in particular, and works at a company whose main product is implemented in FePy.

Just to make it clear, this is not a book for someone who wants to pick up programming from scratch. There's little hand holding here. The obligatory introduction to Python is quickly and efficiently dealt with before moving on to the main subject - accessing .NET classes and creating a GUI application. For those coming from C# to IronPython, there is a whole chapter on unit testing with Python, as well as another that covers everything from Python magic methods to metaprogramming. A particularly nice feature of the authors' style is to logically link each section to the following, so that you always understand the point of what you've just read and where it fits into the bigger picture.

From my point of view, it's nice to see that explanations are provided for those using the free version of Visual Studio, Visual Studio Express, as this means that I can test out all of the examples myself. Other chapters yet to be finalised cover using IronPython in a webbrowser (apparently IronPython can run in Silverlight, the new browser plugin from Microsoft), and extending IronPython with C# or VisualBasic.

In short, for any serious users of IronPython, this book is a must have. It may also convince those using C# to make the switch to a better and more productive life with Python...

Reacting to questions

2008-07-15T13:04:00.000-07:00

Don't blink or you'll miss it. My 15 minutes of fame starts now.* David Bradley tracked me down in cyberspace, and turned me into a Reactive Profile.

*Update: 15 minutes now over.

Pipeline Python - Generate a workflow

2008-07-10T08:49:00.000-07:00

Workflow packages such as Pipeline Pilot, Taverna and KNIME allow the user to graphically create a pipeline to process molecular data. A downside of these packages is that the units of the workflow, the nodes, process data sequentially. That is, no data gets to Node 2 until Node 1 has finished processing all of it. Correction (thanks Egon): The previous line is plain incorrect. Both KNIME and Taverna2, at least, pass on partially processed data as soon as it's available.

Wouldn't it be nicer if they worked more like Unix pipes, that is, as soon as some data comes out of Node 1 it gets passed onto the next Node and so on. This would have three advantages: (1) you get the first result quicker, (2) you don't use up loads of memory storing all of the intermediate results, (3) you can run things in parallel, e.g. Node 2 could start processing the data from Node 1 immediately, perhaps even on a different computer.

Luckily, there is a neat feature in Python called a generator that allows you to create a pipeline that processes data in parallel. Generators are functions that return a sequence of values. However, unlike just returning a list of values, they only calculate and return the next item in the sequence when requested. One reason this is useful is because the sequence of items could be very large, or even infinite in length. (For a more serious introduction, see David Beazley's talk at PyCon'08, which is the inspiration for this blog post.)

Let's create a pipeline for processing an SDF file that has three nodes: (1) a filter node that looks for the word "ZINC00" in the title of the molecule, (2) a filter node for Tanimoto similarity to a target molecule, (3) an output node that returns the molecule title. (The full program is presented at the end of this post.)

# Pipeline Python!
pipeline = createpipeline((titlematches, "ZINC00"),
                            (similarto, targetmol, 0.50),
                            (moltotitle,))

# Create an input source
dataset = pybel.readfile("sdf", inputfile)

# Feed the pipeline    
results = pipeline(dataset)

The variable 'results' is a generator, so nothing actually happens until we request the values returned by the generator...

# Print out each answer as it comes
for title in results:
  print title

The titles of the molecules found will appear on the screen one by one as they are found, just like in a Unix pipe. Note how easy it is to combine nodes into a pipeline.

Here's the full program:


import re
import os
import itertools

# from cinfony import pybel
import pybel

def createpipeline(*filters):
  def pipeline(dataset):
      piped_data = dataset
      for filter in filters:
          piped_data = filter[0](piped_data, *filter[1:])
      return piped_data
  return pipeline

def titlematches(mols, patt):
  p = re.compile(patt)
  return (mol for mol in mols if p.search(mol.title))

def similarto(mols, target, cutoff=0.7):
  target_fp = target.calcfp()
  return (mol for mol in mols if (mol.calcfp() | target_fp) >= cutoff)

def moltotitle(mols):
  return (mol.title for mol in mols)

if __name__ == "__main__":
  inputfile = os.path.join("..", "face-off", "timing", "3_p0.0.sdf")
  dataset = pybel.readfile("sdf", inputfile)
  findtargetmol = createpipeline((titlematches, "ZINC00002647"),)
  targetmol = findtargetmol(dataset).next()

  # Pipeline Python!
  pipeline = createpipeline((titlematches, "ZINC00"),
                            (similarto, targetmol, 0.50),
                            (moltotitle,))

  # Create an input source
  dataset = pybel.readfile("sdf", inputfile)

  # Feed the pipeline    
  results = pipeline(dataset)

  # Print out each answer as it comes through the pipeline    
  for title in results:
      print title

So, if in future someone tells you that Python generators can be used to make a workflow, don't say "I never node that".

Image: Travis S.

Chemoinformatics p0wned by cheminformatics

2008-07-08T01:11:00.000-07:00

The headline says it all. After two weeks, and thousands, er...40 votes, the results are in: 26 for cheminformatics, 14 for that other one.

Next week, bioinformatics versus binformatics...

Cheminformatics toolkit face-off: Speed (Python vs Java vs C++)

2008-07-06T12:20:00.000-07:00

It's a bit meaningless to compares speeds of different toolkits performing the same operations. Functionality is really the reason you are going to favour one toolkit over another. Here I'm focusing on comparing the speed of accessing the same toolkit from CPython or Jython versus accessing it directly in its original language. In other words, what price do you pay to be able to work in Python rather than C++ or Java? (I'll discuss the advantages of working in Python in a later post.)

To begin with, let's consider accessing the CDK from Python versus from Java. The test input file for all of these tests is the first subset of the drug-like ligands in ZINC, 3_p0.0.sdf, which contains 24098 molecules. The three test cases are (1) Iterate over all of the molecules, (2) iterate over all of molecules and write out the molecular weight, (3) calculate 25 descriptor values for the first 228 molecules. I implemented these in Java, and using cinfony. For example, here's the cinfony version for test case 2:

import time
from cinfony import cdk

t = time.time()
for mol in cdk.readfile("sdf", "3_p0.0.sdf"):
   print mol.molwt
print time.time() - t

Here are the results (times are in seconds):

Method	Test 1	Test 2	Test 3
Java	22.2	38.9	31.7
CPython (cinfony, cdkjpype)	34.0	72.6	38.2
Jython (cinfony, cdkjython)	23.7	44.4	34.0

It's clear that accessing the CDK from Jython is almost as fast as using it in a Java program. However, there is an overhead associated with using it from CPython except where, as in Test 3, most of the time is spent in computation.

Next, let's look at accessing OpenBabel from Python versus from C++. Here I will compare the following tests cases: (1) iterate over all of the molecules, (2) iterate over all of molecules and write out the molecular weight, (3) apply 30 steps of a forcefield optimisation to the first 200 molecules. Here's an example of the cinfony script for (2). Notice any similarities to the one above? :-)

import time
from cinfony import pybel

t = time.time()
for mol in pybel.readfile("sdf", "3_p0.0.sdf"):
   print mol.molwt
print time.time() - t

Here are the results (note: measured on a different machine than the tests above):

Method	Test 1	Test 2	Test 3
C++	77.7	126.0	56.8
CPython (cinfony, Pybel)	78.3	132.9	60.0
Jython (cinfony, Jybel)	80.4	135.7	59.4

In short, the cost of using Pybel or Jybel is small.

Technical notes: For the CDK, I calculated the values of all descriptors (except IP) that didn't return an array of values. This came to 25 descriptors. I also skipped over one molecule, #20, that took several seconds to process. Jython can natively access Java libraries such as the CDK, but to access the CDK from CPython, cinfony uses JPype. cinfony uses the Python SWIG wrappers to access OpenBabel from CPython; Jython is using the Java SWIG wrappers for OpenBabel. I need to repeat the runs a few times on a quiet machine to get better figures, but I note that the figures do not include the cost of loading the OpenBabel DLL or starting up the JVM.

Image credit:MacRonin47

ANN: OpenBabel 2.2.0 Released

2008-07-05T01:39:00.000-07:00

Here is the official announcement from Geoff Hutchison:

I am very happy to finally announce the release of Open Babel 2.2.0, the latest stable version of the open source chemistry toolbox.

This release represents a major update and should be a stable upgrade, strongly recommended for all users of Open Babel. Highlights include improved force fields and coordinate generation, conformer searching, enhanced plugins including molecular descriptors, filters, and command-line transformations. Many formats are improved or added, including CIF, mmCIF, Gaussian cube, PQR, OpenDX cubes, and more. Improved developer API and scripting support and many, many bug fixes are also included.

What's new? See the full release notes.

To download, see our Install Page.

For more information, see the project website.

I would like to personally thank a few people for making this release a great one. In alphabetical order, Jean Bréfort, Andrew Dalke, Marcus Hanwell, Chris Morley, Noel O'Boyle, Kevin Shepherd, Tim Vandermeersch, and Ugo Varetto.

This is a community project and we couldn't have made this release without you. Many thanks to all the contributors to Open Babel including those of you who submitted feedback, bug reports, and code.

Cheers,
-Geoff

ANN: cinfony 0.2 - the "easy to install" version

2008-06-25T00:18:00.000-07:00

A new version of cinfony is now available. The good news is that it now much easier to install for Windows users.

All you need now is to have Python and Java, download the CDK and the RDKit, edit a configuration file, and away you go using OpenBabel, the RDKit and the CDK from Python. The full instructions are here.

I've updated some of the docstrings, but documentation is still sparse and the Pybel documentation is still the best guide, as linked to on the cinfony web site.

Here's an example of what's possible with cinfony. Suppose you want to convert a SMILES string to 3D coordinates with OpenBabel, create a 2D depiction of that molecule with the RDKit, calculate descriptors with the CDK, and write out an SDF file containing the descriptor values and the 3D coordinates.

from cinfony import rdkit, cdk, pybel
mol = pybel.readstring("smi", "CCC=O")
mol.make3D()
rdkit.Molecule(mol).draw(show=False,
                         filename="aldehyde.png")
descs = cdk.Molecule(mol).calcdesc()
mol.data.update(descs)
mol.write("sdf", filename="aldehyde.sdf")

This new version of cinfony also includes Jybel, so that you can now access both the CDK and OpenBabel from Jython, as well as from CPython.

The Forced Authorship Licence - Get your users to write papers for you

2008-06-24T06:22:00.000-07:00

You are probably familiar with commercial software licences. You may even have heard of open source licences. But are you familiar with the Forced Authorship Licence (FAL) model? Let me give you an example from real life (the name has been removed to focus on the actual licence):

"X is available free of charge for researchers belonging to Academic community. The download and use of X is subject to the X Academic Licence...Every users is associated with one of the X team. This will help the user in installing and using X and will be co-author of the first paper published by the user, containing X results. You must contact one of the three group leaders and discuss your proposed project before applying for the use of X."

I think this is a great idea. Think of all the publications. If I had thought of this a few years ago, I would now have 30 extra papers instead of the 30 citations of GaussSum.

But it's never to late to start. And why stop at software? If I'm going to be competing with people that use the FAL, I need to think smarter. From now on, all of my papers, software, blog posts, personal communications, and any ideas that arose while reading my papers, attending my talks or reading my posters will carry the Noel O'Blog Licence (NOABL - 'A' for apostrophe and don't you forget it). Instead of citing me, any resulting publication must carry my name as sole author, in bold - no, make that in fire in letters thirty feet high - and when applying for grants, I'm allowed to list these publications together with my own work.

Sounds fair, doesn't it? Oh - I almost forgot. It's spelt "Noel M. O'Boyle". And don't forget that apostrophe.

Image: Martin Deutsch

O No - It's cheminformatics!

2008-06-23T14:56:00.000-07:00

It is time to cast your vote in the greatest polarising debate of our times. Yes indeed: should it be cheminformatics or chemoinformatics? Vote now (see poll on right).

Not to sway any undecided voters, but I'm definitely in favour of "cheminformatics". My main reason is that I'm worried that if the other camp win out, they'll probably decide to change more words: we'll end up doing chemoistry, like our Australian cousins.

You have 13 days to cast your vote...

An RSS feed for the CCL list

2008-06-23T13:44:00.000-07:00

Apparently some people still read the CCL (Computational Chemistry List) using email. I gave up on that some time ago. Here's an RSS feed I threw together some time ago, and which you might find useful: CCL feed

If you don't know what an RSS feed is, and how it might be useful, it's quicker just to test it out than to explain. First of all, you need a feed reader: for example get an account on Google Reader, a free online RSS reader. Next, subscribe to whatever RSS feeds you want, by clicking on "Add subscription" and copying and pasting the URL of the RSS feed into the box that appears.

Here's how the feed is created:

import sys
import email
import datetime
import pdb

from ftplib import FTP
from StringIO import StringIO

import PyRSS2Gen

def breakdaily(messages):
    """
    >>> import pickle
    >>> a = pickle.load(open("tmp20070105"))
    >>> len(breakdaily(a))
    1
    """
    broken = []
    message = []
    for line in messages.split("\n"):
        if line.startswith("From owner-chemistry@ccl.net"):
            broken.append("\n".join(message))
            message = []
        else:
            message.append(line)
    broken.append("\n".join(message))
    return broken[1:]

def getlatest(N):
    ftp = FTP('ftp.ccl.net')
    ftp.login('anonymous')
    ftp.cwd('/pub/chemistry/archived-messages')
    # ftp.dir()

    listoffiles = ftp.nlst()

    # Get current year
    year = datetime.datetime.now().year

    # Get N most recent messages
    messagetot = 0
    months = [str(x).zfill(2) for x in range(1, 13)]
    months.reverse()
    days = [str(x).zfill(2) for x in range(1, 32)]
    days.reverse()
    msgs = []


    while messagetot<N:
        ftp.cwd(str(year))

        for month in months:
            ftp.cwd(month)
            availabledays = ftp.nlst()
            for day in [x for x in days if x in availabledays]:
                # Go thru in reverse order but exclude days that are
                # non-existent
                messages = StringIO()
                ftp.retrbinary("RETR %s" % day, messages.write)
##                pickle.dump(messages.getvalue(), open("tmp%i%s%s" % (year,month,day), "w"))
                listmsgs = breakdaily(messages.getvalue())
                # print "="*24 + "\n", messages.getvalue()
                for i,msg in enumerate(listmsgs):
                    msg_content = email.message_from_string(msg)
                    text = ""
                    for part in msg_content.walk():
                        if (part.get_content_maintype()=="text" and
                            part.get_content_type()=="text/plain"):
                            text = part.get_payload()
                            text = text.replace("\n","<br/>").decode("iso-8859-1", "strict")
                    msgs.append( (year, month, day, msg_content, i+1, text) )
                messagetot += len(listmsgs)
                if messagetot>=N:
                    break
            if messagetot>=N:
                break
            ftp.cwd("..")
        ftp.cwd("..")
        year -= 1 # Continue into the previous year

    ftp.quit()
    msgs.reverse()

    return msgs

def main():

    print "\nStarting..."

    messages = getlatest(100)
##    outputfile = open("messages.pickle", "w")
##    pickle.dump(messages, outputfile)
##    outputfile.close()
##    import sys
##    sys.exit(1)
##    messages = pickle.load(open("messages.pickle", "r"))

    rssitems = []
    for year, month, day, msg, id, messagetext in messages:

        # Add the new item
        newitem = PyRSS2Gen.RSSItem(
                 title = msg['Subject'],
                 link = "http://ccl.net/cgi-bin/ccl/message-new?%d+%s+%s+%s" % (
                         year, month, day, str(id).zfill(3)),
                 description = messagetext,
        ## What's a guid? A globally unique id...used by RSS readers
        ## to determine whether they've seen a particular news item
        ## already
                 guid = PyRSS2Gen.Guid("http://ccl.net/cgi-bin/ccl/message-new?%d+%s+%s+%s" % (
                         year, month, day, str(id).zfill(3))),
                 pubDate = msg['Date']
                 )
        rssitems.append(newitem)

    rss = PyRSS2Gen.RSS2(
        title = "CCL",
        link = "http://www.ccl.net",
        description = "RSS Feed of the world's greatest computational chemistry "
                      "mailing list, chemistry@ccl.net (the CCL list)",
        lastBuildDate = datetime.datetime.now(),
        items = rssitems)
    rss.write_xml(open("ccl.rss", "w"))

    print "Finishing...\n"

def test():
    import doctest
    doctest.testmod()

def do_debugger(type, value, tb):
    pdb.pm()

if __name__=="__main__":

    # sys.excepthook = do_debugger
    main()

    ## test()

Jmol gets competition - OpenAstexViewer now available

2008-06-16T23:19:00.000-07:00

AstexViewer has just gone open source (LGPL), and is available from http://openastexviewer.net/.

Both Jmol and AstexViewer are 3D chemical structure applets that run in the browser. AstexViewer was originally developed by Mike Hartshorn at Astex Therapeutics for in-house visualisation of protein crystal structures. Although available at no cost for some time, it has not until now been open source.

Jmol is in many ways an open source success story. It has several enthusiastic developers (who make new releases with new features every few days!), a very busy mailing list, a large number of users worldwide, and even a recent book. It will be interesting to see whether OpenAstexViewer is sufficiently different to attract users away from this well-established project.

In any case, here's to diversity. Hopefully both projects will get interesting ideas from each other.

Cheminformatics toolkit face-off - SMARTS matching

2008-05-21T13:23:00.000-07:00

SMARTS strings are really useful for substructure searching in molecules. They are sort of like regular expressions for molecules. The idea and syntax of SMARTS comes from the people at Daylight.

Rajarshi Guha, together with Dazhi Jiao (Uni. Indiana), have put together a test suite for SMARTS matchers, and run the test suite against the CDK, OpenBabel and OEChem. Greg contributed results for RDKit. The overall performance is described on Rajarshi's website.

Here's the current summary of the results for the 158 test cases, but check Rajarshi's page for a more up-to-date picture:

Toolkit	True	False
CDK	149	9
OpenBabel	149	9
RDKit	151	7
OpenEye	150	8

Image: Chris Radcliff

Cheminformatics toolkit face-off - Depiction Part 2

2008-05-20T12:31:00.000-07:00

One of the aims of the previous toolkit face-off was to get some feedback on the best options for drawing images with different toolkits. I also hoped that a comparison of the images would allow bugs to be easily identified. And finally, I thought that a bit of competition might help improve depiction and structure diagram generators across the board.

Since the last post:

Molinspiration have approached me to be included in the face-off
I've learnt that I should remove hydrogens from CDK and OASA depictions and control the size of the generated image
Beda has done amazing work enhancing depiction in OASA, and has in fact now released OASA separately from BKChem as an independent cheminformatics toolkit
Greg is just about to release a new version of the RDKit which has improved depiction, and he has fixed the bug identified by PMR
And I've fixed some bugs myself in cinfony

Now the images themselves (same dataset as before): [depiction] and [structure diagram generation].

Notes:

The face-off involves the open source toolkits the CDK, OASA and the RDKit, and the proprietary toolkits Cactvs and molinspiration.
If you want to help improve these depictions or coordinate generators, why not leave a comment below suggesting specific ways to improve them, or highlighting specific things they could do better.
Double bond stereochemistry doesn't seem to be preserved by the CDK, but this is possibly my fault (I'm awaiting a reply to an email to the cdk-users mailing list)
Several people complained to me last time that I didn't give OASA a fair chance. In order to make it up, this time round I've made OASA the star attraction. The coordinates generated by the three open source toolkits are all depicted by OASA

RSS feeds for chemistry projects on SourceForge

2008-05-12T00:00:00.000-07:00

Wouldn't it be nice to know whether any of your favourite chemistry projects has released a new version? Or to keep abreast of the latest registrations of chemistry projects on SourceForge? No? I guess it's just me then.

Anyway, here are some RSS feeds I've thrown together to allow me to do just that:

latest chemistry releases
registrations of chemistry projects: latest and two months old

Note: the two-months-old RSS feed is better if you want to find some actual code or a working website when you click through.

Want to do the same for another software category? Here's the code (requires BeautifulSoup and PyRSS2Gen):


import datetime
import urllib

from BeautifulSoup import BeautifulSoup
import PyRSS2Gen

def download():
    urls = ["http://sourceforge.net/search/index.php?words=trove%3A%28384%29"
            "&sort=latest_file_date&sortdir=desc&offset=0&limit=100&"
            "type_of_search=soft&pmode=0",
            "http://sourceforge.net/search/index.php?words=trove%3A%28384%29"
            "&sort=registration_date&sortdir=desc&offset=0&limit=100&"
            "type_of_search=soft&pmode=0"]
    urllib.urlretrieve(urls[0], "releases.html")
    urllib.urlretrieve(urls[1], "registrations.html")

def converttodate(text):
    if text=="(none)":
        return None
    t = map(int, text.split("-"))

    return datetime.datetime(*t)

def makerss(filename, items, sortby, title):
    rss = PyRSS2Gen.RSS2(
        title = title,
        link = "http://baoilleach.blogspot.com/2008/05/rss-feeds-for-chemistry-projects-on.html",
        description = "baoilleach's RSS feed of "
                      "Chemistry projects on SourceForge",

        lastBuildDate = datetime.datetime.now(),

        items = [
           PyRSS2Gen.RSSItem(
             title = item["title"],
             link = item["link"],
             description = item["description"],
             guid = PyRSS2Gen.Guid("%s %s" % (item["title"], item['lastrelease'])),
             pubDate = item[sortby])
           for item in items]
           )

    rss.write_xml(open(filename, "w"))

def analyse(project):
    ans = {}
    ans['title'] = project.a.string
    ans['link'] = "http://sf.net" + project.a['href']

    data = project.parent.parent
    ans['lastrelease']  = converttodate(data('td')[5].string.strip())

    ans['registered'] = converttodate(data('td')[4].string.strip())

    data = project.parent.parent.findNextSibling()
    ans['description'] = data.td.contents[0].strip()
    if not ans['description']:
        ans['description'] = data.td.contents[2].strip()

    return ans

def processfile(filename):
    html = open(filename, "r").read()
    soup = BeautifulSoup(html)
    projects = soup.findAll(lambda tag: tag.name=="h3" and tag.a
                            and tag.a['href'].startswith("/projects/"))

    data = [analyse(project) for project in projects]
    return data

if __name__=="__main__":
    download()

    data = processfile("registrations.html")
    sometimeago = datetime.datetime.now() - datetime.timedelta(days=60)
    olddata = [d for d in data if d['registered'] <= sometimeago]
    makerss("oldregistrations.rss", olddata, "registered", "Registrations 60 days ago on SF")
    makerss("newregistrations.rss", data, "registered", "Latest registrations on SF")

    data = processfile("releases.html")
    makerss("latestreleases.rss", data, "lastrelease", "Latest releases on SF")

Review of "Gnuplot in Action"

2008-05-08T07:04:00.001-07:00

"Famous scientific plotting package". This succinct yet accurate description of gnuplot appears on the gnuplot SF project page. Despite being one of the most well-known open source scientific projects around, it has taken 15 years for the first book on gnuplot to be published. "Gnuplot in Action" by Phillip K. Janert will be published by Manning Publications later this year (one chapter available free) and aims to give an overview of how to use gnuplot, describing everything from how to read in data, style graphs, save images to more advanced topics such as macros, using gnuplot for CGI or calling it from scripting languages.

In the past I have used gnuplot to get a quick look at data in text files, and also to look at surfaces. I am used to looking up the gnuplot help system to figure things out, but although it's exhaustive, it's also a bit exhausting especially when you are simply trying to figure out whether something is possible or not. The fact that the gnuplot site has a long list of frequently-asked questions is also a sign of poor documentation. So there's definitely room for a book of this type.

Although I was afraid this book would simply be a long list of commands and switches (of which gnuplot has many, to be sure), "Gnuplot in Action" takes a much more interesting approach, illustrating everything with examples and focusing on the most useful options while simply referring the reader to the gnuplot help for more obscure options. The author isn't afraid of sympathising with the user, e.g. "Nevertheless, it is very different from the behavior we have come to expect from user interfaces in most programs" (in reference to how images are saved in gnuplot). Although the author has been using gnuplot for 15 years, he obviously is still aware of what the most annoying features are, and in this case, describes a handy macro to look after the horrible mess of writing an image to disk and resetting the terminal afterwards.

Somewhat bizarrely, the blurb for the book is pitched at business users (analysts and "Six Sigma Black Belts") and doesn't refer to scientists at all. But don't worry, the examples in the book contain everything from share price analysis to looking at phase transitions (reflecting the author's background in theoretical physics as well as in business), and in any case, the nature of the examples doesn't affect the tone of the book in any way.

So, overall I liked this book. I particularly liked the sidebars and comments on particular techniques. For example, the section on 3d plots includes a short discussion of the pitfalls of using 3d plots and alternative plotting methods that might be more suitable. The section on color, includes a similar warning as well a discussion on palette design (which references Why Should Engineers and Scientists Be Worried About Color?). This turns the book from being a straightforward technical manual to something more of a discussion of techniques. There are actually a couple of chapters still to come (13, 14, 15) which appear to be focused on where/how to use particular graphical analysis techniques.

Any negatives? Well, I wasn't very interested in the section on reading data from disk and transforming it - I think that anyone who knows Perl/Python is going to reshape/filter their data before getting it into gnuplot, but admittedly the book claims to be suitable for non-programmers. Also, the author's favorite interactive terminal, wxt, is not available for Windows and so I was a bit disappointed when I fired up gnuplot and couldn't find it. Regarding the connection between gnuplot and Python, the author doesn't mention the widely-used Gnuplot.py (although I'm not a big fan of this).

But overall, these are minor quibbles. I recommend this book to anyone who wants to make the most of gnuplot, that "famous scientific plotting package".

Cheminformatics toolkit face-off - Depiction

2008-05-07T23:50:00.000-07:00

Now for some pretty pictures as well as some not so pretty. Yes, it's the turn of the structure diagram generators (SDGs) to strut their stuff and throw some shapes. How do they perform for 100 random compounds from PubChem?

Here are my results for depiction and structure diagram generation (Note: I will move these links to my regular hosting in the near future). The images were generated using cinfony, the code is here and the dataset is here.

Some comments:
(0) Rich Apodaca has written an overview of Open Source SDGs.
(1) 2D coordinate generation is independent of depiction. A SDG typically has both parts but coordinates could be generated with one toolkit and depicted with another.
(2) Looking good is not the same as chemical accuracy. But looking good is important too! :-)
(3) OASA (Obviously Another Stupid Acronym) is a Python SDG which is part of BKchem by Beda Kosata. To depict images, it uses Pycairo. OASA currently does not support stereochemistry (e.g. across double bonds) and so the generated coordinates will not be chemically-accurate in some circumstances. This of course does not affect its ability to depict coordinates generated by other toolkits. OASA can depict wedges and bonds but I haven't yet implemented this in cinfony.
(4) It is not possible to use the CDK's depiction mechanism natively from CPython using JPype (it's fine from Jython though). That is why I use OASA to depict the CDK's generated coordinates. Technically, this is because JPype doesn't allow subclassing of Java classes (JPanel in this case). It does however, allow Python classes to implement a Java interface. It would be great if the CDK could add a convenience interface to allow this (I can supply more details off-blog).
(5) The PubChem images and coordinates are generated by the Cactvs toolkit.
(6) There are two sets of RDKit images. Since the current release has quite poor depictions, I decided to also include the results from a development branch, which uses Aggdraw.
(7) It is important to consider how to handle hydrogens. With OASA, I just drew all the hydrogens. This is probably not a good idea.
(8) It is likely that I have not used the best parameters, etc. for generating some of these images. I welcome any suggestions to improve image quality.
(9) OASA wasn't able to layout CID 1373132, and so that molecule was not included (I guess I should have just caught the exception and continued). The error message was:

Traceback (most recent call last):
   File "cdkjpype.py", line 464, in draw
     oasa.cairo_out.cairo_out().mol_to_cairo(mol, filename)
   File "bkchem\oasa\oasa\cairo_out.py", line 85, i
n mol_to_cairo
     self._draw_edge( e)
   File "bkchem\oasa\oasa\cairo_out.py", line 121,
in _draw_edge
     side += reduce( operator.add, [geometry.on_which_
side_is_point( 
start+end, (
self.transformer.transform_xy( a.x,a.y))) for a
  in ring if a!=v1 and a!=v2])
   File "bkchem\oasa\oasa\geometry.py", line 92, in
  on_which_side_is_point
     b = atan2( y2-y1, x2-x1)
ValueError: math domain error

(10) PubChem entries with more than 1 connected component were not included in this test. (As a result, the number of molecules shown is actually less than 100.)

Image: ARTS

ANN: cinfony 0.1 - the toolkit of cheminformatics toolkits

2008-05-07T11:12:00.001-07:00

What?
cinfony presents a common Python API to three open source cheminformatics toolkits: OpenBabel, RDKit and the CDK. cinfony makes it easy to carry common cheminformatics tasks, but allows you to access the underlying toolkit objects at any point.

Why?
(1) In the past, cheminformaticians have chosen a particular toolkit and have only used that toolkit. This is because of the impedance mismatch between APIs of different toolkits, and indeed, the language differences between different toolkits (CDK is Java, OpenBabel and RDKit are C++, for example). cinfony makes it easy to start using a new toolkit because the API for reading/writing files, calculating fingerprints, creating 2D depictions, etc. is the same no matter whether you are using OpenBabel, RDKit or the CDK.

(2) Different toolkits have different capabilities. cinfony allows you to access all of these capabilities in the same Python script. For example, you can read a molecule using OpenBabel, draw it using RDKit and calculate descriptors using the CDK.

Where?
cinfony 0.1 is now available for download from googlecode.

Who?
Early adopters. This is pretty much a test release.

Why are there so many dependencies?
Oops, look at the time. Gotta go. Have fun installing.

Image: DryIcons (modified)

Shortened SMILES for URLS - Why a big SMILES need not mean a long face

2008-04-30T18:59:00.000-07:00

Rajarshi has put together a RESTful web service that returns a 2D depiction of a molecule (using the CDK) given a SMILES string. However, it turns out that the SMILES has a maximum size of 255 characters (appears to be an OS or mod_python limit). That should be enough, sez you, but for a random sample of 100 PubChem molecules, 9 are greater than the limit (these include explicit hydrogens, I should point out).

Between Rajarshi and myself, we worked out a solution: zip the SMILES and then base64 it (so that it can be used in a URL). Here are the lengths of the original SMILES strings, the bzipped2 strings, and finally the base64ed strings:

280 101 136
316 113 152
432 110 148
282 106 144
320 92 124
372 115 156
452 143 192
326 96 128
268 94 128

In each case, the final string is less than half the size of the original string. The actual strings themselves all begin with the same few letters (don't ask me why) which means that the same RESTful URL can be used to handle both the SMILES strings and their compressed cousins. For example, sildenafil can be viewed at both this url (containing the SMILES) and this url (containing the encoded form).

Notes: (1) Python methods bz2.compress() and base64.urlsafe_b64encode() are used. (2) Bzip2 was used instead of gzip or zip simply because it has an easier to use API. (3) A SMILES string needs to have a certain length before the procedure will actually result in any compression. In the sildenafil example, the encoded string is in fact longer than the original (108 vs. 61)

Dear LazyWeb - How to write a unittest for when a Python module is unavailable?

2008-04-30T05:12:00.000-07:00

I have a Python module (cinfony.pybel) which I am testing with the standard unittest module. One of the classes in the module has an optional dependency; if the Python Imaging Library (PIL) is available, a particular draw() method will work, otherwise it should raise an ImportError with a helpful message (I am simplifying the situation somewhat but this is the general problem).

I have two related questions:

(1) Should I place the "import Image" statement for PIL at the start of the draw() method, or at the start of the module (and catch the exception but set a flag)? Personally, I feel that its more Pythonic to keep all of the imports at the start of the module; that is, to be up-front about dependencies, optional or not.

(2) How can I test this code in a unittest? That is, how can I call this draw() method with PIL available, and a second time, without PIL available?

Cheminformatics toolkit face-off - LogP and TPSA

2008-04-23T19:46:00.000-07:00

Updated 23/04/08 (see below)

In a recent post, I discussed some thoughts I had about using the Pybel API as a generic API for cheminformatics toolkits. I've started to implement this (more on this at a later date) and have begun writing some tests to compare results from different toolkits. To begin with, I was hoping that RDKit (Jan2008), CDK (1.0.2) and OpenBabel (dev, soon-to-be 2.2.0) could agree on LogP and TPSA values calculated with the same group contribution approach.

The code is below, but first here are the results:

ALogP
=====
****Example1*****
The calculated value from JCICS, 1999, 39, 868 is 1.400000
cinfony.rdk {'MolLogP': 1.4007999999999998,
             'pyMolLogP': 1.4007999999999998}
cinfony.pybel {'LogP': 0.80749999999999988}



****Example2*****
The calculated value from JCICS, 1999, 39, 868 is 2.750000
cinfony.rdk {'MolLogP': 2.7486000000000006,
             'pyMolLogP': 2.7486000000000006}
cinfony.pybel {'LogP': 1.6415999999999995}



****methane*****
The calculated value from JCICS, 1999, 39, 868 is 0.144100
cinfony.rdk {'MolLogP': 0.6331, 'pyMolLogP': 0.6331}
cinfony.pybel {'LogP': 0.14410000000000001}



TPSA
====
****metoprolol*****
The calculated value from JMC, 2000, 43, 3714 is 50.700000
cinfony.cdk {'tpsa': 50.719999999999999}
cinfony.rdk {'TPSA': 50.719999999999999}
cinfony.pybel {'TPSA': 50.719999999999999}



****nordiazempam*****
The calculated value from JMC, 2000, 43, 3714 is 41.500000
cinfony.cdk {'tpsa': 41.460000000000001}
cinfony.rdk {'TPSA': 41.460000000000001}
cinfony.pybel {'TPSA': 41.460000000000001}

The good news is that they all agree on the TPSA of the two examples. The bad news is that CDK 1.0.2 does not calculate ALogP, and OpenBabel doesn't seem to be doing well for the more complicated examples (which are taken from the JCICS paper). It's possible that OpenBabel is using a different parameterisation than the paper I reference. I'm not sure either if there's any difference between the two ALogP values that RDKit calculates. Hopefully, the people that know will comment below.

Update (23/04/08): (1) I calculated the LogP of methane incorrectly - it should be 0.63610. (2) RDKit has a transposition error in the value for C in methane (only a minor error). (3) OpenBabel calculates the correct answer if hydrogens are added. I'm going to recommend that the descriptor adds the hydrogens itself, to avoid user errors in future. (Done! Thanks to Chris Morley.) (4) The two ALogP descriptors in RDKit are implementations in two different languages, one in C++, the other in Python. They should be identical. (Thanks to Greg Landrum of RDKit for identifying several of these issues)

Here's the code...


from cinfony import cdk, rdk, pybel


descs = [("ALogP", "JCICS, 1999, 39, 868"),
         ("TPSA", "JMC, 2000, 43, 3714")]
toolkits = [(cdk, {'ALogP': None, 'TPSA': ['tpsa']}),
            (rdk, {'ALogP': ['MolLogP', 'pyMolLogP'],
             'TPSA': ['TPSA']}),
            (pybel, {'ALogP': ['LogP'], 'TPSA': ['TPSA']})]
testcases = {'ALogP':
                {"methane": (0.1441, 'C'),
                 "Example1": (1.40, 'c1(O)ccccc1(OC)'),
                 "Example2": (2.75, 'c1ccccc1c1ccccn1')},
             'TPSA':
                {"metoprolol": (50.7,
             "c1(OCC(O)CNC(C)C)ccc(CCOC)cc1"),
                 "nordiazempam": (41.5,
                 "c1c(Cl)ccc2NC(=O)CN=C(c3ccccc3)c12")}}

for desc, ref in descs:
    print desc + "\n" + "="*len(desc)
    for k, v in testcases[desc].iteritems():
        print "****%s*****" % k
        print "The calculated value from %s is %f" % (ref, v[0])
        for toolkit in toolkits:
            if toolkit[1][desc]:
                mol = toolkit[0].readstring("smi", v[1])
                print toolkit[0].__name__,
                print mol.calcdesc(toolkit[1][desc])
        print "\n\n"

Image credit: pfly

Cheminformatics toolkit face-off - Molecular weight

2008-04-22T08:19:00.000-07:00

Next up in this series of toolkit comparisons is the molecular weight. Here, cinfony is going to support two ideas:

molwt: the standard molar mass given by IUPAC atomic masses (amu). This is calculated using OpenBabel's OBMol::GetMolWt(), CDK's MFAnalyser.getCanonicalMass, and the MolWt descriptor in RDKit.
exactmass: the mass given by the most abundant isotope. This is available through OBMol::GetExactMass(), CDK's MFAnalyser.getMass, and is not available in RDKit.

The gold standard is the element data in the Blue Obelisk Data Repository (BODR), which was created for just this purpose.

The following are the results for a hydrogen atom and a carbon atom according to the BODR and each of the toolkits. The first figure in parentheses is the molwt, then the exactmass.

BODR: (1.00794, 1.007823032)
     (12.0107, 12)

Pybel: (1.00794, 1.007825032)
      (12.0107, 12.0)

RDKit: (1.008)
    (12.011)

CDK: (1.0079400539398193, 2.0156500339508057)
    (12.010700225830078, 16.031299591064453)

OpenBabel is the only toolkit that gets everything right. RDKit is doing okay, but should consider using BODR data in future. The CDK presents the most intriguing results. I'm not sure whether the Python/Java interface has introduced noise into the molwt values, but they are exactly in agreement with the BODR up to the 7th decimal place or so, after which something goes weird. On the other hand, it's quite clear that the CDK's exactmass is simply not behaving as advertised. It is using Deuterium for the hydrogen and C-16 for the carbon. (I note that MFAnalyser has already been replaced in the CDK development code.)

Image credit: bugmonkey

Noel O'Blog celebrates one year at ACS New Orleans

2008-04-03T01:32:00.001-07:00

To celebrate the first anniversary of Noel O'Blog, I'm going to be attending the ACS National Meeting in New Orleans, so say hi if you see me. I'll mostly be hanging around the CINF and COMP sessions, or at the CCDC (Cambridge Crystallographic Data Centre) stand in the Exposition Hall, and for sure, at the CINF reception on Tuesday evening.

If you're tired of hearing talks about how docking doesn't work, etc., etc., come along early on Monday and I'll explain one way to improve a scoring function (in the docking software GOLD) to pick out those pesky actives. Any other ideas people might have would be very welcome.

Pybel as a generic API for cheminformatics libraries - proof of concept using CDK

2008-03-28T10:40:00.000-07:00

I'm very interested in interoperability of open source chemistry codes. Following a comment by Egon on a recent post of mine, I started wondering whether the Pybel API could used with other cheminformatics libraries as a backend.

The advantage of this for the user would be (a) to reduce the learning curve - if you know how to use Pybel, you can access any of several different cheminformatics libraries with the same syntax, (b) the same scripts could be used to carry out a particular analysis using different cheminformatics libraries - different libraries may have different fingerprints, descriptors or implementations of particular algorithms (this is of course also useful for cross-checking the results of different programs) and (c) help reduce the divide between different cheminformatics toolkits (interoperability!!).

The rationale behind Pybel (described in the paper) lends itself to this use. Pybel doesn't attempt to wrap all the functionality of OpenBabel, but only the most common tasks in cheminformatics. For advanced options, or additional functionality, you can go behind the scenes and access OpenBabel directly. As a result, I propose that the Pybel API represents a generic API (one of many possible, of course) for accessing any cheminformatics library.

To test this, I have created CDKabel, a proof of concept which shows that the Chemistry Development Kit (CDK) can be accessed using Pybel syntax through Jython. CDKabel does not yet pass all of the Pybel tests, but there's enough to show that the approach has some merit. Compare the following: here's some Python code using Pybel and OpenBabel:

C:\Documents and Settings\oboyle>python25
Python 2.5 (r25:51908, Sep 19 2006, 09:52:17) [MSC v.1310 32
bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more inf
ormation.
>>> from pybel import *
>>> for mol in readfile("sdf", "head.sdf"):
...         print "Molecule has molwt of %.2f and %d atoms" %
                  (mol.molwt, len(mol.atoms))
...
Molecule has molwt of 122.12 and 15 atoms
Molecule has molwt of 332.49 and 28 atoms
>>>

Now here's some Jython code with CDKabel and CDK:

D:\Tools\CDK>set CLASSPATH=cdk-1.0.2.jar

D:\Tools\CDK>..\jython2.2.1\jython
Jython 2.2.1 on java1.6.0_05
Type "copyright", "credits" or "license" for more informa
tion.
>>> from cdkabel import *
>>> for mol in readfile("sdf", "head.sdf"):
...     print "Molecule has molwt of %.2f and %d atoms" %
               (mol.molwt, len(mol.atoms))
...
Molecule has molwt of 122.04 and 15 atoms
Molecule has molwt of 331.96 and 28 atoms
>>>

Well, at least they agree on the number of atoms :-) (It's my fault - CDK has like, ten different ways of calculating the molecular mass, and I just chose randomly :-) )

I've only spent a few minutes throwing CDKabel together, so it doesn't do much beyond the example shown. However, if interested, you can download it and try it for yourself.

I'd appreciate comments on the idea that there is a core Python API that could be usefully applied to several cheminformatics libraries. Would anyone use CDKabel if it were available?