Get 3D distance matrix with rdkit #RDKit

I updated rdkit of my env from 20180301 to 20180303 with anaconda. ;-)
When I want to get 3D distance matrix of the molecule I use Get3DDistanceMatrix method.
But I found that rdDistGeom.GetMoleculeBoundsMatrix returns almost same results.
3DDistance matrix is useful for feature of 3D QSAR.
I would like to use these method.

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import rdDistGeom as molDG

mol = Chem.MolFromSmiles('CCC')
bm = molDG.GetMoleculeBoundsMatrix(mol)
bm
Out[]: 
array([[0.        , 1.524     , 2.51279063],
       [1.504     , 0.        , 1.524     ],
       [2.43279063, 1.504     , 0.        ]])

AllChem.EmbedMolecule(mol)
Out[]: 0

dm=AllChem.Get3DDistanceMatrix(mol)
dm
Out[]: 
array([[0.        , 1.50401361, 2.47848679],
       [1.50401361, 0.        , 1.50913087],
       [2.47848679, 1.50913087, 0.        ]])

mol = Chem.MolFromSmiles('c1ccncc1')

bm = molDG.GetMoleculeBoundsMatrix(mol)
bm
Out[]: 
array([[0.        , 1.38925641, 2.42894217, 2.80158084, 2.42894217,
        1.38925641],
       [1.36925641, 0.        , 1.38925641, 2.39940019, 2.81851281,
        2.42894217],
       [2.34894217, 1.36925641, 0.        , 1.35507278, 2.3697344 ,
        2.81851281],
       [2.68158084, 2.31940019, 1.33507278, 0.        , 1.35507278,
        2.39940019],
       [2.34894217, 2.64739923, 2.2897344 , 1.33507278, 0.        ,
        1.38925641],
       [1.36925641, 2.34894217, 2.64739923, 2.31940019, 1.36925641,
        0.        ]])

IAllChem.EmbedMolecule(mol)
Out[31]: 0

dm=AllChem.Get3DDistanceMatrix(mol)
dm
Out[]: 
array([[0.        , 1.38579572, 2.3376015 , 2.70976216, 2.34726355,
        1.38818637],
       [1.38579572, 0.        , 1.39124965, 2.30040952, 2.68623726,
        2.37414223],
       [2.3376015 , 1.39124965, 0.        , 1.35847264, 2.28377475,
        2.61995998],
       [2.70976216, 2.30040952, 1.35847264, 0.        , 1.35652537,
        2.3746108 ],
       [2.34726355, 2.68623726, 2.28377475, 1.35652537, 0.        ,
        1.3905889 ],
       [1.38818637, 2.37414223, 2.61995998, 2.3746108 , 1.3905889 ,
        0.        ]])

And also I would like to extract 3D pharmacophore feature from molecule.
Example is below.

import os
from rdkit import Geometry
from rdkit import RDConfig
from rdkit.Chem import AllChem
from rdkit.Chem import ChemicalFeatures
from rdkit.Chem.Pharm3D import Pharmacophore
FEAT = os.path.join(RDConfig.RDDataDir, "BaseFeatures.fdef")
featfact = ChemicalFeatures.BuildFeatureFactory(FEAT)
mol = Chem.MolFromSmiles('c1cccnc1')
AllChem.EmbedMolecule(mol)
feats = featfact.GetFeauresFromMol(mol)
for feat in feats:
    ...:     print(feat.GetFamily())
    ...:     pos = feat.GetPos()
    ...:     print(pos.x, pos.y, pos.z)
    ...:     
    ...:     
Acceptor
-1.042491325161116 -0.7212337219142626 -0.7060687598321927
Aromatic
2.7755575615628914e-16 -1.942890293094024e-16 -0.2252839133689341

RDKit has many useful features for chemoinformatics!

Published by iwatobipen

I'm medicinal chemist in mid size of pharmaceutical company. I love chemoinfo, cording, organic synthesis, my family.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.

%d bloggers like this: