Pybel - Hack that SD file

At some point, most cheminformaticians need to hack at SD files. Maybe they need to extract the values of particular data fields, or add in some new ones. Unfortunately, our friends 'awk', 'grep' and 'Excel' (!) are not very useful for this so we need to use a scripting language.

The examples shown here use the Python module Pybel, included with the Open Babel distribution. Although it's a Python module, all of the hard work is done by the underlying C++ library.

Convert an SD file to a spreadsheet

Given an SD file with calculated descriptor values in the fields, I want to write out a Tab-separated file containing a SMILES string in the first column, the title in the second column, and then the descriptor values in the remaining columns.

import pybel

inputfile = pybel.readfile("sdf", "ace_ligands.sdf")
#  We need to read the first molecule to find
#  the descriptor names
mol = inputfile.next()
print "\t".join(["SMILES", "Title"] + mol.data.keys())

# Print out information on the first molecule
print "\t".join([mol.write("smi").split()[0], mol.title]
                 + mol.data.values())
# Print out information on the remaining molecules
for mol in inputfile:
    print "\t".join([mol.write("smi").split()[0], mol.title]
                 + mol.data.values())

Notes:

• The data fields in an SD file are in the 'data' attribute of a molecule, which behaves like a dictionary with keys() and values().
• The "\t".join(mylist) concatenates a list of strings around Tab characters.
• If you convert a molecule to a SMILES string, the result is a string containing the SMILES string, a Tab, the molecule title, and a carraige return. To simply access the SMILES string, you need to split() it and access the first element.

Add Rule-of-5 descriptor values to an SD file

Although Lipinksi's Rule-of-5 (or fives) may be both a Sacred Cow and an Evil Metric (according to the Great Molecular Crapshoot), we still love it. Here's how to add RO5 descriptor values to an SD file (requires Open Babel 2.1.1):

import pybel

HBD = pybel.Smarts("[!#6;!H0]")
HBA = pybel.Smarts("[$([$([#8,#16]);!$(*=N~O);" +
                   "!$(*~N=O);X1,X2]),$([#7;v3;" +
                   "!$([nH]);!$(*(-a)-a)])]")

def moredesc(mol):
    ans = {}
    ans['RotBonds'] = mol.OBMol.NumRotors()
    ans['HBD'] = len(HBD.findall(mol))
    ans['HBA'] = len(HBA.findall(mol))
    ans['molwt'] = mol.molwt
    return ans

output = pybel.Outputfile("sdf", "LipinskiRulesOK.sdf")

for mol in pybel.readfile("sdf", "NoRules.sdf"):
    mol.data.update(mol.calcdesc())
    mol.data.update(moredesc(mol))
    output.write(mol)

output.close()

Notes:

• The 'calcdesc' method of a Molecule returns a dictionary containing descriptor values for LogP, Molecular Refractivity (MR), and Polar Surface Area (PSA).
• The 'update' method of a dictionary adds the contents of one dictionary to another.