[Rdkit-discuss] Pharmacophore matching

[Hilleke M. <mat...@ph...>]

Dear all,

I am new to working with the Pharm3D module and I have a question regarding pharmacophore matching: Is it possible to match a molecule to a pharmacophore that requires one atom to match to two different features (e.g. H-bond acceptor and donor)? Here is a dummy code example for a case where I am trying to match ethanol back to all of its pharmacophore features (I am using rdkit version 2022.09.4):


import os
from rdkit import Chem, RDConfig
from rdkit.Chem import AllChem, ChemicalFeatures, rdDistGeom
from rdkit.Chem.Pharm3D import Pharmacophore, EmbedLib

feature_definition_file = os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef')
featFactory = ChemicalFeatures.BuildFeatureFactory(feature_definition_file)

m1 = Chem.MolFromSmiles('CCO') # ethanol
m1 = Chem.AddHs(m1)
AllChem.EmbedMolecule(m1)

features = featFactory.GetFeaturesForMol(m1)

# This is to keep track of features
for i in features:
    print(i.GetFamily())
#features = features[1:]    #removing donor feature will make it work
print()
for i in features:
    print(i.GetFamily())

pcophore = Pharmacophore.Pharmacophore(features)
m2 = Chem.MolFromSmiles('CCO') # ethanol again
m2 = Chem.AddHs(m2)
AllChem.EmbedMolecule(m2)

bool_match, matches = EmbedLib.MatchPharmacophoreToMol(m2, featFactory, pcophore)
print("bool_match =", bool_match)

boundsMat = rdDistGeom.GetMoleculeBoundsMatrix(m2)
failed, boundsMatMatched, matched, matchDetails = EmbedLib.MatchPharmacophore(matches, boundsMat, pcophore, useDownsampling=True)
print("matched =", matched)



MatchPharmacophore() appears to be working only by removing one of the two features derived from the oxygen. Am I doing something wrong or is there a neat way to fix this?

Kind regards,
Mattis