Pybel - Just how unique are your molecules? Part II

I recently posted on using SMILES, Fingerprints and InChIs to identify unique molecules in a dataset. Geoff pointed out that I should have been using canonical SMILES ("can" in Open Babel), instead of non-canonical SMILES ("smi" in Open Babel). Why?

Well, the same molecule will always have the same canonical SMILES. Whereas depending on the order of the atoms in the input file, the same molecule might have a different non-canonical SMILES. Since I was interested in identifying unique molecules, I should have been using canonical SMILES.

Method

The same as before, just with "can" instead of "smi". The code used is only slightly different from that in the previous post.

Results

Here is the output of the script (see code below):

Starting...

The number of molecules is 12712

Are there any molecules with the same FP, SMI and InChI?

There are 2 molecules with 2 duplicates
There is 1 molecule with 162 duplicates
There is 1 molecule with 661 duplicates
There is 1 molecule with 1098 duplicates


The number of (unique) molecules is 10792

Are there any remaining molecules with the same canonical 
   SMILES?

There are 815 molecules with 2 duplicates
There are 61 molecules with 3 duplicates
There are 36 molecules with 4 duplicates
There are 10 molecules with 5 duplicates
There are 9 molecules with 6 duplicates
There are 9 molecules with 7 duplicates
There is 1 molecule with 8 duplicates
There are 3 molecules with 10 duplicates
There are 3 molecules with 11 duplicates
There are 2 molecules with 13 duplicates
There is 1 molecule with 14 duplicates
There is 1 molecule with 24 duplicates
There is 1 molecule with 34 duplicates

Are there any remaining molecules with the same fingerprint?

None found

Are there any remaining molecules with the same InChI?

There are 2 molecules with 2 duplicates


...Finished!

Discussion

If you compare the results here with those in my previous post, you will find that that there are two instances where two molecules have the same FP, InChI and non-canonical SMILES, but they have different canonical SMILES! By eye the two molecules are identical. This appears to be a bug in the canonicalisation routine and I've reported it to OpenBabel. The details are:

non-canonical SMILES is CC1OC(C[NH](C1)CC(C(F)(F)F)O)C but canonical SMILES is...

• ZINC03883383: C[C@H]1O[C@@H](C)C[NH](C1)C[C@H](O)C(F)(F)F
• ZINC03883386: C[C@@H]1O[C@H](C)C[NH](C1)C[C@H](O)C(F)(F)F

non-canonical SMILES is C1C(CC(C(C1O)O)O)(C(=O)O)O but canonical SMILES is...

• ZINC03870192: O[C@@H]1CC(O)(C[C@H](O)C1O)C(=O)O
• ZINC03870194: O[C@H]1CC(O)(C[C@@H](O)C1O)C(=O)O

Although my original aim was simply to give an example of how to use Pybel, this trivial analysis of ZINC has identified areas for improvement in both ZINC and OpenBabel. Without access to such a large dataset as ZINC, errors such as those identified here would go unnoticed. If anyone has any more ideas for tests, let me know...

Pybel - Just how unique are your molecules?

Cheminformaticians often want to remove or identify duplicates in a dataset of molecules. There are a couple of popular ways to do this: ensure that each molecule has a unique SMILES string, has a unique Daylight-type fingerprint, or has a unique InChI.

But what's the difference? Is there any? Does a unique fingerprint imply a unique InChI, and vice versa?

Method

Enough with the questions. Let's just get some numbers. Our test set is the first slice of the drug-like subset of the ZINC database. Specifically, it's the file 3_p0.0.mol2.gz. According to ZINC, this should contain "only a single form of each molecule". We're going to use Pybel, the Python module included with Open Babel, to calculate the three common measures of uniqueness and check for duplicates.

Results

Here is the output of the script (see code below):

Starting...

The total number of molecules is 12712

Are there any molecules with the same FP, SMI and InChI?
There are 4 molecules with 2 duplicates
There is 1 molecule with 162 duplicates
There is 1 molecule with 661 duplicates
There is 1 molecule with 1098 duplicates

The total number of (unique) molecules is 10790

Are there any remaining molecules which have the same FP \
    and SMI?

There are 846 molecules with 2 duplicates
There are 55 molecules with 3 duplicates
There are 23 molecules with 4 duplicates
There are 3 molecules with 5 duplicates
There is 1 molecule with 7 duplicates

Are there any remaining molecules which have the same InChI?

...Finished!

Discussion

So it seems that uniqueness is not quite as easy to measure as you might think. Even once you remove all molecules that share the same fingerprint, SMILES string, and InChI, there are still many that share the same fingerprint and SMILES string, although the InChIs are unique.

The results shown here either expose errors in OpenBabel, or possible errors in the creation of the ZINC database. In either case, this is an example of two Open chemistry resources being used to test and validate each other. Hopefully these results will aid in the development of even better software or data.

Code

import os
import glob
import pybel
import cPickle

def duplicates(keys, data):
    """Find the duplicates in 'data' in terms of particular
    keys.

    'keys' is a list containing some selection of the
    indices 1, 2 and 3.
    """
    ans = {}
    for d in data:
        ans.setdefault(tuple([d[key] for key in keys]),[]) \
            .append(d[TITLE])
    result = {}
    for k, v in ans.iteritems():
        if len(v) > 1: # len(v) is the number of duplicates
            result[len(v)] = result.get(len(v), 0) + 1
    for k in sorted(result.keys()):
        if result[k]>1:
            print "There are %d molecules with %d" \
                  " duplicates" % (result[k], k)
        else:
            print "There is %d molecule with %d" \
                  " duplicates" % (result[k], k)
    return ans

def readfrommol2(filename):
    data = []
    for mol in pybel.readfile("mol2", filename):
        data.append(( mol.title,
            "%s" % mol.calcfp('FP2').bits,
            mol.write("smi").split()[0],
            mol.write("inchi").rstrip() ))
    outputfile = open("dups.pickle", "w")
    cPickle.dump(data, outputfile)
    outputfile.close()
    return data

def readfrompickle():
    inputfile = open("dups.pickle", "r")
    data = cPickle.load(inputfile)
    inputfile.close()
    return data

if __name__ == "__main__":
    print "Starting...\n"

    if not os.path.isfile("dups.pickle"):
        data = readfrommol2("3_p0.0.mol2")
    else:
        data = readfrompickle()

    TITLE, FP, SMI, INC = range(4)

    # How many molecules are there?
    print "The number of molecules is %d\n" \
            % len(data)

    # Find molecules with the same FP, SMI and INC
    print "Are there any molecules with the same" \
          " FP, SMI and InChI?"
    answer = duplicates( [SMI, FP, INC], data)

    # Select only unique molecules
    uniquedata = [(v[0], k[0], k[1], k[2]) for k, v \
            in answer.iteritems()]
    print "\nThe number of (unique) molecules is %d\n" \
            % len(uniquedata)

    # Find molecules with the same FP and SMI but not INC
    print "Are there any remaining molecules which have" \
          " the same FP and SMI?\n"
    duplicates( [SMI, FP], uniquedata )

    # Find any remaining molecules with the same InChI
    print "\nAre there any remaining molecules which have" \
          " the same InChI?\n"
    duplicates( [INC], uniquedata ) # Nothing is found

    print "...Finished!"