For each record from the benzodiazepine file, print the total number of heavy atoms in each record (that is, exclude hydrogens). You must get the data from the connection table and not from the tag data. (In particular, don't use PUBCHEM_HEAVY_ATOM_COUNT .)
The output is one output line per record, containing the count as an integer. If at all possible, show how to read directly from the gzip'ed input SD file. The first 5 heavy atom counts are 21, 24, 22, 22 and 22. A copy of the expected output is available for reference.
Some of the hydrogens in the data set have atomic weights of 2 or 3, which may cause difficulties for toolkits which convert 1H hydrogens to implicit form while leaving 2H and 3H to explicit form. Please make sure your output agrees with the reference output.
CDK/Java
Save the code to "countHeavyAtoms.java" then do
javac countHeavyAtoms.java java -cp ./:$HOME/cdk-1.2.4.1.jar countHeavyAtoms > countHeavyAtoms.out
import java.io.*;
import org.openscience.cdk.interfaces.IChemObjectBuilder;
import org.openscience.cdk.io.iterator.IteratingMDLReader;
import org.openscience.cdk.Molecule;
import org.openscience.cdk.DefaultChemObjectBuilder;
import org.openscience.cdk.interfaces.IAtom;
public class countHeavyAtoms {
public static void main(String[] args) {
FileReader sdfile = null;
try {
/* CDK does not automatically understand gzipped files */
sdfile = new FileReader(new File("benzodiazepine.sdf"));
} catch (FileNotFoundException e) {
System.err.println("benzodiazepine.sdf not found");
System.exit(1);
}
IteratingMDLReader mdliter = new IteratingMDLReader(sdfile,
DefaultChemObjectBuilder.getInstance());
Molecule mol = null;
while (mdliter.hasNext()) {
mol = (Molecule) mdliter.next();
int numHeavies = 0;
for (IAtom a : mol.atoms()) {
if (a.getAtomicNumber() > 1) {
numHeavies += 1;
}
}
System.out.println(numHeavies);
}
}
}
CDK/Cinfony/Jython
from cinfony import cdk
manip = cdk.cdk.tools.manipulator.AtomContainerManipulator()
for mol in cdk.readfile("sdf", "benzodiazepine.sdf"):
print len(manip.getHeavyAtoms(mol.Molecule))
Indigo
#include "base_cpp/scanner.h"
#include "base_cpp/output.h"
#include "molecule/molecule.h"
#include "molecule/sdf_loader.h"
#include "molecule/molfile_loader.h"
int main (int argc, char *argv[])
{
const char *idfield = 0;
bool allow_query = false;
int i;
if (argc <= 1)
{
fprintf(stderr,
"Usage: count_heavy file.sdf [-id idfield]\n");
return -1;
}
for (i = 2; i < argc; i++)
{
if (strcmp(argv[i], "-id") == 0)
{
if (++i == argc)
{
fprintf(stderr, "expecting an identifier after '-id'\n");
return -1;
}
idfield = argv[i];
}
else if (strcmp(argv[i], "-allow-query") == 0)
allow_query = true;
else
{
fprintf(stderr, "unknown parameter: %s\n", argv[i]);
return -1;
}
}
const char *filename = argv[1];
try
{
FileScanner scanner(filename);
SdfLoader loader(scanner);
int cnt = 0;
if (idfield != 0)
loader.initProperties(idfield);
while (!loader.isEOF())
{
loader.readNext();
const char *id = 0;
if (idfield != 0 && loader.properties.at2(idfield) != 0)
id = loader.properties.at(idfield).ptr();
cnt++;
if (idfield == 0)
printf("molecule #%d: ", cnt);
else if (id == 0)
printf("molecule #%d, no %s: ", cnt, idfield);
else
printf("molecule #%d, %s=%s: ", cnt, idfield, id);
fflush(stdout);
Molecule mol;
try
{
BufferScanner molscan(loader.data);
MolfileLoader molload(molscan);
molload.loadMolecule(mol, allow_query);
int heavy_count = 0;
for (i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
if (mol.getAtom(i).label > ELEM_H)
heavy_count++;
printf("%d heavy atoms\n", heavy_count);
}
catch (Exception &e)
{
printf("error: %s\n", e.message());
continue;
}
}
}
catch (Exception &e)
{
fprintf(stderr, "error: %s\n", e.message());
return -1;
}
return 0;
}
Instructions:
- Unpack 'graph' and 'molecule' projects into some folder
- Create 'utils' folder nearby
- Paste the following code into utils/cansmiles.cpp file:
- Compile the file using the following commands:
$ cd graph; make CONF=Release32; cd ..
$ cd molecule; make CONF=Release32; cd ..
$ cd utils
$ gcc count_heavy.cpp -o count_heavy -O3 -m32 -I.. -I../common ../molecule/dist/Release32/GNU-Linux-x86/libmolecule.a ../graph/dist/Release32/GNU-Linux-x86/libgraph.a -lpthread -lstdc++ - Run the program like that:
$ ./count_heavy benzodiazepine.sdf -id pubchem_compound_cid
OpenBabel/Python
import openbabel as ob
obconversion = ob.OBConversion()
obconversion.SetInFormat("sdf")
obmol = ob.OBMol()
notatend = obconversion.ReadFile(obmol, "benzodiazepine.sdf.gz")
while notatend:
print obmol.NumHvyAtoms()
# reuse the molecule for 0.3 seconds of savings vs. making an new one
obmol.Clear()
notatend = obconversion.Read(obmol)
OpenBabel/Pybel
This builds on and simplifies OpenBabel
import pybel
for mol in pybel.readfile("sdf", "benzodiazepine.sdf.gz"):
print mol.OBMol.NumHvyAtoms()
OpenEye/Python
from openeye.oechem import *
ifs = oemolistream()
ifs.open("benzodiazepine.sdf.gz")
for mol in ifs.GetOEGraphMols():
print OECount(mol, OEIsHeavy())
# Can't use
# print mol.NumAtoms()
# because the input structure contains explicit hydrogens.
RDKit/Python
I could also have iterated over the atoms and only counted the non-hydrogen ones, but I doubt that would be any faster or elegant.
from rdkit import Chem
# RDKit does not seem to support gzip'ed SD files
suppl = Chem.SDMolSupplier("benzodiazepine.sdf")
heavy_patt = Chem.MolFromSmarts("[!#1]")
for mol in suppl:
print len(mol.GetSubstructMatches(heavy_patt))
The method "mol.GetNumAtoms()" is supposed to report only the non-hydrogen counts, but actually reports the number of explicit atoms in the molecular graph.
RDKit/Cinfony/Python
from cinfony import rdk
for mol in rdk.readfile("sdf", "benzodiazepine.sdf"):
print mol.calcdesc(['HeavyAtomCount'])['HeavyAtomCount']
Unix command-line
It is possible to do this on the Unix command-line using a single (very long) line. To make it easier to understand I've written it using an awk script named 'heavy_atom_counts.awk'. Use it as:
gzcat benzodiazepine.sdf.gz | awk -f heavy_atom_counts.awk
# This file is named "heavy_atom_counts.awk"
# Check track of the line number in the record
{lineno += 1}
# Get the atom counts from the record
lineno == 4 { num_atoms = int(substr($0, 1, 3)); num_hydrogens = 0}
# Count the number of hydrogens
4 < lineno && lineno <= 4+num_atoms && substr($0, 32, 2) == "H " {num_hydrogens++}
# And the end of the record, print the number of heavy atoms and reset
/^\$\$\$\$/ {print num_atoms - num_hydrogens; lineno = 0}