Re: [Rdkit-discuss] pharmacophore align

[Greg L. <gre...@gm...>]

Dear Taka,

On Thu, Jul 19, 2012 at 9:44 AM, Taka Seri <ser...@gm...> wrote:
> Dear all.
> I have some questions about Pharmacophore.
> I want to align molecules that from SDF files,  by Pharmacophore.
> At first, I set Pharmacophore by using "EmbedPharmacophore()" method.
> And checked molecules by "EmbedPharmacophore()" method.
> Then,  generated 3D structures by "EmbedPharmacophore()" method.
> Finally, I want to align these conformers to the Pharmacophore.
> Could anyone give me a small example?

There is unfortunately no really good sample code available for this.
There is some very old code from a GUI application that was built with
the RDKit but that is no longer supported (or part of the current SVN
app) here:
http://rdkit.svn.sourceforge.net/viewvc/rdkit/trunk/Python/qtGui/Search3D/SearchUtils.py?revision=2&view=markup&pathrev=5
in the function AlignMatchToReference()


Here's an attempt to distill that information down:

First you need the alignment package:

from rdkit.Numerics  import rdAlignment

which contains the function GetAlignmentTransform:

In [5]: rdAlignment.GetAlignmentTransform?
Type:       function
Base Class: <type 'builtin_function_or_method'>
String Form:<Boost.Python.function object at 0x102fe0dc0>
Namespace:  Interactive
Docstring:
GetAlignmentTransform( (object)refPoints, (object)probePoints [,
(object)weights=[] [, (bool)reflect=False [, (int)maxIterations=50]]])
-> object :
    Compute the optimal alignment (minimum RMSD) between two set of points


     ARGUMENTS:

        - refPoints : reference points sepcified as a N by 3 Numeric array or
                      sequence of 3-sequences or sequence of Point3Ds
        - probePoints : probe points to align to reference points - same format
                      restrictions as reference points apply here
        - weights : optional numeric vector or list of weights to
associate to each pair of points
        - reflect : reflect the probe points before attempting alignment
        - maxIteration : maximum number of iterations to try to minimize RMSD

     RETURNS:

        a 2-tuple:
          - SSD value for the alignment
          - the 4x4 transform matrix, as a Numeric array


This will give the transformation required to align one set of points
(the ph4 points from your embedded molecule) to another set of points
(the ph4 points from your reference ph4).

To do this you need the positions of each of the features:
probePts = [list(x.GetPos()) for x in probeFeats]
refPts = [list(x.GetPos()) for x in refFeats]

And then you align them:
ssd,tform = Aligner.GetAlignmentTransform(refArr,probeArr,weights=weights)

If your molecule has no chiral centers, it's a good idea to try
reflecting the alignment:
ssd2,tform2 = Aligner.GetAlignmentTransform(refArr,probeArr,weights=weights,reflect=True)
if ssd2<ssd:
 	tform = tform2
 	ssd = ssd2


the two return values are the sum of squared deviations of the atomic
positions (easily convertible to the RMSD) and the transformation
matrix. You can apply the transformation matrix to your embedded probe
molecule (to align it to the ph4) with the function
AllChem.TransformMol:
AllChem.TransformMol(probeMol,tform)

I hope this helps.

-greg