Heavy atom counts from an SD file

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.

Suggestion: Please add a couple of '$$$$' SD data lines to the file, to trip those naive implementations which think that a '$$$$' is automatically a record separator!

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

CDK/Groovy
Save the code to "countHeavyAtoms.groovy" then do:

groovy countHeavyAtoms.groovy > countHeavyAtoms.out

Indigo/C++
Instructions: $ 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++ $ ./count_heavy
 * 1) Unpack 'graph' and 'molecule' projects into some folder
 * 2) Create 'utils' folder nearby
 * 3) Paste the above code into utils/count_heavy.cpp file
 * 4) Compile the file using the following commands:
 * 1) Run the program like that:

OpenBabel/Pybel
This builds on and simplifies OpenBabel

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.

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.

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

{lineno += 1}
 * 1) This file is named "heavy_atom_counts.awk"
 * 2) Check track of the line number in the record

lineno == 4 { num_atoms = int(substr($0, 1, 3)); num_hydrogens = 0}
 * 1) Get the atom counts from the record

4 < lineno && lineno <= 4+num_atoms && substr($0, 32, 2) == "H " {num_hydrogens++}
 * 1) Count the number of hydrogens

/^\$\$\$\$/ {print num_atoms - num_hydrogens; lineno = 0}
 * 1) And the end of the record, print the number of heavy atoms and reset

Cactvs/Tcl
molfile loop "benzodiazepine.sdf.gz" eh { puts [ens get $eh E_HEAVY_ATOM_COUNT] }

The toolkit can read compressed files directly. The loop construct automatically discards the structure read for each record, so this works on arbitrarily big files. Also, the format of the input file can be anything which can be identified automatically, so this is not limited to SD files.

Cactvs/Python
for eh in Molfile('benzodiazepine.sdf.gz'): print(eh.E_HEAVY_ATOM_COUNT)

The Python language version has the same capabilities as the Tcl variant.