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!"