Visualize pharmacophore with RDKit #RDKit #Pymol #Chemoinformatics

Some years ago, I wrote a post about how to communicate pymol and RDKit. In the post, I demonstrated how to visualize Phamarcophore in rdkit.

And recently I got a query about the post and think about more efficient way.
What I want to write in the post is new approach to visualize pharmacophore with RDKit.

To do it, I ran pymol in server mode. The environment of pymol session is not need to same in rdkit env. I installed pymol via homebrew. So this version is python 3.7. And my environment of rdkit is python3.6.

iwatobipen$ pymol -R

Then write code on jupyter notebook. I picked up several mols from cdk2.sdf which is provided from rdkit Docs/Book/data folder.

from IPython import display
import os
import subprocess
from rdkit import Chem
from rdkit import RDConfig
from rdkit.Chem import AllChem
from rdkit.Chem import rdMolAlign
from rdkit.Chem import rdShapeHelpers
from rdkit.Chem import rdMolDescriptors
from rdkit.Chem import PyMol
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
import copy
import pprint
mols = Chem.SDMolSupplier('./before_align_cdk2.sdf', removeHs=False)
cdk2mol = [m for m in mols]
for m in cdk2mol:
    AllChem.EmbedMolecule(m, AllChem.ETKDGv2())

Then aligned molecules with crippenO3A method. I used first molecule as reference. And after align I saved molecule to align_cdk2.sdf file.

cdk2mol2 = copy.deepcopy(cdk2mol)
crippen_contribs = [rdMolDescriptors._CalcCrippenContribs(mol) for mol in cdk2mol2]
ref = cdk2mol2[0]
ref_contrib = crippen_contribs[0]
targets = cdk2mol2[1:]
targets_contrib = crippen_contribs[1:]

for i, target in enumerate(targets):
    crippenO3A = rdMolAlign.GetCrippenO3A(target, ref, targets_contrib[i], ref_contrib)
    crippenO3A.Align()

w = Chem.SDWriter('./align_cdk2.sdf')
w.write(ref)
for mol in targets:
    w.write(mol)
w.close()

Key function of ShowFeatures is ‘rdkit/Chem/Features/ShowFeats.py’
To use the function, it is needed to run the code like below.

iwatobipen$ python ShowFeats.py 

I want to call the function in jupyter notebook. So I used subprocess.
Following example shows before alignment and after alignment.

showfeatpath = os.path.join(RDConfig.RDCodeDir, 'Chem/Features/ShowFeats.py')

At first visualize molecules before align. ‘display’ function can display png image on notebook.

# Before align
v = PyMol.MolViewer()
v.DeleteAll()
process = subprocess.Popen(['python', showfeatpath, '--writeFeats','./before_align_cdk2.sdf'], stdout=subprocess.PIPE)
stdout = process.communicate()[0]
png=v.GetPNG()
display.display(png)

Hmm, the function can detect pharmacophores of each molecule but not aligned form. Next visualize molecules

#After align
v = PyMol.MolViewer()
v.DeleteAll()
process = subprocess.Popen(['python', showfeatpath,'--writeFeats','./align_cdk2.sdf'], stdout=subprocess.PIPE)
stdout = process.communicate()[0]
png=v.GetPNG()
display.display(png)

Now I could get following image.

And I set –writeFeats option. The option provides position of each pharmacophores.

res = stdout.decode('utf-8').replace('\t', ' ').split('\n')
pprint.pprint(res)

['# Family    X     Y    Z    Radius  # Atom_ids',
 'Donor -4.447 -0.580 -2.127 1.0 # 11',
 'Donor -2.369 -0.730 -2.637 1.0 # 13',
 'Donor -4.121 -0.009 0.275 1.0 # 14',
 'Donor -3.602 0.536 2.585 1.0 # 17',
 'Acceptor 2.288 -0.534 -1.549 1.0 # 5',
 'Acceptor -0.161 -0.294 -0.684 1.0 # 7',
 'Acceptor -2.369 -0.730 -2.637 1.0 # 13',
 'Acceptor -4.121 -0.009 0.275 1.0 # 14',
 'Acceptor -1.919 0.112 0.908 1.0 # 16',
 'PosIonizable -3.325 -0.576 -2.046 1.0 # 9 13 12 11 10',
 'Aromatic -3.325 -0.576 -2.046 1.0 # 9 10 11 12 13',
 'Aromatic -2.826 -0.102 -0.038 1.0 # 8 9 10 14 15 16',
 'Hydrophobe 3.473 -0.083 0.406 1.0 # 2',
 'LumpedHydrophobe 3.899 0.272 0.319 1.0 # 2 1 3',
 'Donor -4.491 -0.607 -2.118 1.0 # 2',
 'Donor -2.390 -0.682 -2.613 1.0 # 5',
 'Donor -4.131 -0.045 0.283 1.0 # 9',
 'Donor -3.541 0.496 2.576 1.0 # 10',
 'Acceptor -2.390 -0.682 -2.613 1.0 # 5',
 'Acceptor -1.879 0.136 0.883 1.0 # 7',
 'Acceptor -4.131 -0.045 0.283 1.0 # 9',
 'Acceptor -0.163 -0.211 -0.758 1.0 # 11',
 'Acceptor 2.410 -1.330 -1.015 1.0 # 17',
 'PosIonizable -3.355 -0.566 -2.032 1.0 # 3 2 1 5 4',
 'Aromatic -3.355 -0.566 -2.032 1.0 # 1 2 3 4 5',
 'Aromatic -2.828 -0.097 -0.047 1.0 # 3 4 6 7 8 9',
 'Donor -4.513 -0.557 -2.079 1.0 # 2',
 'Donor -2.426 -0.726 -2.623 1.0 # 5',
 'Donor -4.148 -0.000 0.288 1.0 # 9',
 'Donor -3.550 0.529 2.554 1.0 # 10',
 'Donor 2.549 0.697 -1.275 1.0 # 18',
 'Acceptor -2.426 -0.726 -2.623 1.0 # 5',
 'Acceptor -1.881 0.097 0.872 1.0 # 7',
 'Acceptor -4.148 -0.000 0.288 1.0 # 9',
 'Acceptor -0.137 -0.319 -0.735 1.0 # 11',
 'Acceptor 4.371 2.116 -1.804 1.0 # 17',
 'PosIonizable -3.375 -0.565 -2.021 1.0 # 3 2 1 5 4',
 'Aromatic -3.375 -0.565 -2.021 1.0 # 1 2 3 4 5',
 'Aromatic -2.824 -0.106 -0.053 1.0 # 3 4 6 7 8 9',
 '']

I want to do with Features.ShowFeats function directly from ipython notebook but it was difficult. So I call this function from other process.

BTW, it is interesting and cool code for visualize pharmacophores. RDKit is nice tool for chemoinformatics.

All code is uploaded to my repo.
https://nbviewer.jupyter.org/github/iwatobipen/playground/blob/master/viz_p4core/display_pharmacophore.ipynb
https://github.com/iwatobipen/playground/tree/master/viz_p4core