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!

Advertisements