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

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

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

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: Atomium by bugmonkey (CC BY-NC 2.0)

Noel O'Blog celebrates one year at ACS New Orleans

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

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?

MadMol - Chemistry-aware code

There are definitely better things you can do with the functional group fingerprint (FP4) of OpenBabel, but instead, I've created MadMol. Basically, give MadMol a SMILES string such as CC(=O)CCl and it will tell you all you ever wanted to know about that molecule:

That there's a Primary_carbon. Isn't that a Alkylchloride? Wake up and smell the Ketone, cause that's what it is. Most people wouldn't realise this is a C_ONS_bond. It's a 1,3-Tautomerizable (or your money back). I don't believe it, it's a Rotatable_bond! Wake up and smell the CH-acidic, cause that's what it is.

And here's the code (requires Pybel, and the file SMARTS_InteLigand.txt):

import sys
import random

import pybel

def readsmartsfile(filename="SMARTS_InteLigand.txt"):
   patterns = []
   inputfile = open(filename, "r")
   for line in inputfile:
       line = line.strip()
       if line and line[0]!="#":
           colon = line.find(":")
           name = line[:colon]
           smarts = line[colon+1:].strip()
           patterns.append([pybel.Smarts(smarts), name])
   return patterns

phrases = ["I don't believe it, it's a %s!",
   "Isn't that a %s?",
   "It's a whadyamacallit, a %s.",
   "Looks like a %s to me.",
   "That there's a %s.",
   "Most people wouldn't realise this is a %s.",
   "It's a %s (or your money back).",
   "Wow, a %s. Last time I saw one of these, I hit the"
   "fire alarm and ran.",
   "Could be a %s...yes, I'm sure of it.",
   "It's a %s if I've ever seen one.",
   "Wake up and smell the %s, cause that's what it is.",
   "It's a %s. I wish I had one.",
   "You've hit the jackpot, you and your %s!",
   "A %s. You know, back in the day, we used to have"
   "fun with these.",
   "It takes me back years, this %s does."]

if __name__=="__main__":
   if not len(sys.argv)==2:
       sys.exit("You need a SMILES string like CC(=O)CCl")
   molecule = pybel.readstring("smi", sys.argv[1])
   print "So you want me to tell you about %s, do you?\n" % (
           sys.argv[1],)
   patterns = readsmartsfile()
   print " ".join([random.choice(phrases) % name for
     smarts, name in patterns if smarts.findall(molecule)])

Python - the scripting language of chemistry

Python is a cross-platform scripting language which is easy to learn, and encourages readability and elegant code. If you're a chemist, it's also the most widely-used scripting language out there:

• Cheminformatics
• OpenBabel (Pybel), RDKit, OEChem, Daylight (PyDaylight), Cambios Molecular Toolkit, Frowns, PyBabel and MolKit (both part of MGLTools)
• Computational chemistry
• OpenBabel, cclib, QMForge, GaussSum, PyQuante, NWChem, Maestro/Jaguar, MMTK
• Visualisation
• CCP1GUI, PyMOL, PMV, Zeobuilder, Chimera, VMD, (Avogadro coming soon)
• Scientific programming
• NumPy (interface to ATLAS, LAPACK), Python can interface to C/C++ and FORTRAN, matplotlib (graphing), VTK (3D visualisation, used by CCP1GUI)

Leave a comment below if you feel I've missed something. Here are some dated but related links: ChemPython.org, the text of Andrew Dalke's talk from PyCon 2004, and Python for Chemistry in 21 Minutes.