Dear Fergal,
On Tue, Apr 30, 2013 at 8:51 PM, Fergal Duffy <[email protected]> wrote:
> I would like some help on using the Pharmacophore matching code in RDKit.
Just a standard warning about this: the code in rdkit.Chem.Pharm3D is
under-documented and based on an unpublished algorithm. But then,
you've certainly discovered that already. ;-)
There is a bit of documentation, written as part of the docs for a GUI
tool that no longer works (and is no longer part of the RDKit
distribution) in: $RDBASE/Docs/Programs/RDPharm3D/MethodOverview.htm.
I will, at some point, take a look at pulling this out and integrating
it into the rest of the RDKit docs.
> I am trying to screen molecules using 3D pharmacophores, and I am not
> totally sure I'm going about it correctly. I start with a known
> 3D structure of a ligand, and want to use EmbedMol() to try and find
> molecules that can be made fit the molecule derived from this.
>
> To build a pharmacophore from chemical features in a 3D molecule, I have
> been doing something like:
>
> from rdkit import Chem
> from rdkit import Geometry
> from rdkit import RDConfig
> from rdkit.Chem import AllChem
> from rdkit.Chem.Pharm3D import Pharmacophore, EmbedLib
>
> m1 = Chem.MolFromMolFile('example.sdf')
> FEATURE_DEF_FILE = os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef')
> feat_factory = ChemicalFeatures.BuildFeatureFactory(FEATURE_DEF_FILE)
> feats = feat_factory.GetFeaturesForMol(m1)
> pcophore = Pharmacophore.Pharmacophore(feats)
That's a pharmacophore that includes every feature in the molecule.
This is probably not what you want.
It's more common to include only the features that are (or that you
believe are) important for the interaction you're studying.
>
> I'm reasonably confident that this is correct. Where I am having trouble is
> with what exactly the bounds matrix is supposed to represent,
> and how to supply the atomMatch parameter to EmbedMol.
>
> I am assuming something like:
>
> bounds_matrix = AllChem.GetMoleculeBoundsMatrix(m1)
> m2 = Chem.MolFromSmiles('c1ccccc1C(=O)N')
> can_match, atomMatch = EmbedLib.MatchPharmacophoreToMol(m2, feat_factory,
> pcophore)
> if can_match:
> EmbedLib.EmbedPharmacophore(m2, atomMatch, pcophore, bounds=bounds_matrix)
>
>
> Which should add a conformer to m2. Is this the best way to do it?
Looks right to me. You shouldn't need to include the bounds matrix as
an argument; if you leave it out, the code will generate the matrix
for you. This is is somewhat less efficient if you are embedding the
same molecule multiple times, but it frees you from having to worry
about the details of getting the bounds matrix right.
You should also pay attention to the results of MatchPharmacophore. It
can fail and you want to detect that.
You might want to look into using EmbedLib.EmbedOne() instead of
EmbedPharmacophore(). EmbedOne() does some additional work to clean up
the geometry of the embedded conformations.
>
> Also, what is the correct way to use excluded volumes. I originally assumed
> that they were implemented as a sphere with
> an x,y,z coordinate centre, and a certain radius, but that does not seem to
> be the case from looking at the documentation.
The excluded volume implementation in Pharm3D never did work
particularly well; I don't think I would trust it. Having said that:
yes, an exclusion sphere is defined by a position (relative to the
pharmacophore feature points) and a radius.
I hope this helps at least somewhat,
-greg
------------------------------------------------------------------------------
Introducing AppDynamics Lite, a free troubleshooting tool for Java/.NET
Get 100% visibility into your production application - at no cost.
Code-level diagnostics for performance bottlenecks with <2% overhead
Download for free and get started troubleshooting in minutes.
http://p.sf.net/sfu/appdyn_d2d_ap1
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
