Given a set of drug molecules from Chembl, I need to calculate pairwise distances between the drugs, for clustering. Noel O'Blog helped me do it using RDKit by Greg Landrum, by writing the Python script below. Given Chembl ids. for the input molecules, it calculates pairwise similarities between drugs based on calculating an ECFP4 fingerprint for each drug molecule, and then calculating the distances between the fingerprints for each pair of drugs. The pairwise similarities given in the output file are scaled from 0-1, so can be converted to distances using 1 - similarity. The input file contains duplicates that are distinguished only by undefined/defined stereocentres; these are likely the same molecule, so the script removes duplicates that are distinguished only by stereochemistry.
Thanks Noel!
import collections
from rdkit import Chem, DataStructs
from rdkit.Chem import Descriptors, AllChem
smiles_lookup = dict( (y, x) for (x, y) in (z.split() for z in open("C:\Tools\Topliss\chembl20\chembl_20.ism") if len(z.split())==2))
def readdata(filename):
"""Return dictionary from chemblid to SMILES"""
it = open(filename)
header = next(it)
data = collections.defaultdict(set)
for line in it:
broken = line.rstrip().split()
if broken[2] in data: # Check for duplicates
assert False
data[broken[2]] = broken[1]
return data
def getmol(smi):
"""Keep a single largest component if disconnected"""
mol = Chem.MolFromSmiles(smi)
if "." in smi:
frags = list(Chem.GetMolFrags(mol, asMols=True))
frags.sort(key=lambda x:x.GetNumHeavyAtoms(), reverse=True)
mol = frags[0]
return mol
# The file contains duplicates that are distinguished only by undefined/defined stereocenters
# These are likely the same molecule. Let's remove duplicates that are distinguished only
# by stereochemistry
def mergeStereoisomers(chemblids):
smis = [data[chemblid] for chemblid in chemblids_filtered]
mols = [getmol(smi) for smi in smis]
normalised = {}
for i, mol in enumerate(mols):
smi = Chem.MolToSmiles(mol)
normalised[smi] = chemblids[i] # Stores only one
return normalised.values()
if __name__ == "__main__":
# ECFP setup
fingerprinter = lambda m: AllChem.GetMorganFingerprintAsBitVect(m, 2, nBits=16384)
srcs = ["SMILES4CMPDS_smiles.tsv"]
for src in srcs:
filename = src
print "Looking at %s" % filename
data = readdata(filename)
chemblids = data.keys()
print "Distinct chemblids in original TSV", len(chemblids)
# filter out those molecules without structures in ChEMBL
# ...this step doesn't apply here
chemblids_filtered = chemblids
print "Chembl Ids associated with SMILES", len(chemblids_filtered)
chemblids_final = mergeStereoisomers(chemblids_filtered)
print "Chembl Ids after removing stereoisomers", len(chemblids_final)
smis = [data[chemblid] for chemblid in chemblids_final]
mols = [getmol(smi) for smi in smis]
fps = [fingerprinter(mol) for mol in mols]
N = len(chemblids_final)
with open("%s_pairwise_sim.txt" % src, "w") as f:
# header
f.write("\t".join(chemblids_final) + "\n")
for i in range(N):
row = []
fpA = fps[i]
for j in range(N):
if i==j:
sim = 1.0
else:
fpB = fps[j]
sim = DataStructs.FingerprintSimilarity(fpA, fpB)
row.append("%.2f" % sim)
f.write(chemblids_final[i] + "\t" + "\t".join(row) + "\n")
Find any compounds in a new set ('set 2') that are not very similar to those in my original set ('set 1')
A second problem I had was to see if any compounds that I had found from a PubChem search ('set 2') were very different from a list I had already curated from various sources ('set 1'). Noel O'Blog came to the rescue, helping me to use RDKit to find molecules that are structurally different, defined by calculating similarities of <=0.4, when similarities were calculated using the Tanimoto coefficient of ECFP4 fingerprints as implemented in the RDKit toolkit, http://www.rdkit.org). Here's the Python script, which takes 'set 2' in input file 'pubchem_219_compounds', and 'set 1' in input file 'tmp.csv':
import csv
from rdkit import Chem, DataStructs
from rdkit.Chem import Descriptors, AllChem
def readPCCompounds():
return [x.rstrip().split("\t") for x in open("pubchem_219_compounds")]
def isDigit(char):
return char >= '0' and char <= '9'
def readX20():
ans = []
reader = csv.reader(open("tmp.csv"))
# Skip header
for N in range(7):
reader.next()
for row in reader:
if row[0] == "": break
if row[6][0]=='-': continue
if row[6].startswith("!"):
broken = row[6][1:].split("!")
if len(broken) == 2:
smiles = broken
else:
smiles = [broken[i] for i in range(1, len(broken), 2)]
for x in smiles:
ans.append( [row[0], x] )
else:
ans.append( [row[0], row[6]] )
return ans
def getmol(smi):
"""Keep a single largest component if disconnected"""
mol = Chem.MolFromSmiles(smi)
if "." in smi:
frags = list(Chem.GetMolFrags(mol, asMols=True))
frags.sort(key=lambda x:x.GetNumHeavyAtoms(), reverse=True)
mol = frags[0]
return mol
def getFP(smi):
mol = getmol(smi)
return AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=16384)
if __name__ == "__main__":
pubchems = readPCCompounds()
knowns = readX20()
pubchemfps = [getFP(smi) for name, smi in pubchems]
knownfps = [getFP(smi) for name, smi in knowns]
# for each PubChem entry...
for pubchem, pubchemfp in zip(pubchems, pubchemfps):
print "Considering PubChem entry %s..." % pubchem[0],
# ...compare to each known...
results = []
for known, knownfp in zip(knowns, knownfps):
sim = DataStructs.FingerprintSimilarity(pubchemfp, knownfp)
if sim > 0.4:
results.append( (sim, known[0]) )
if results:
results.sort(reverse=True)
print ", ".join("%s (%.2f)" % (x[1], x[0]) for x in results)
else:
print "nothing found with sim > 0.4"
No comments:
Post a Comment