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
Thank you for this post.
I have come across a similar one here: http://w grep SCORE /mnt/lustre/users/bdiallo/PFDOCK2/Docking/Grim/DrugBank/3vi2_1_protein_DB00152_all* | sort -k 3ww.ub.edu/cbdd/?q=content/mapping-chemical-features-molecules-using-rdkit.
I am under windows and struggling to run Pymol in server mode. I have a WSL and will try to run pymol in server mode from it.
Your post here (https://iwatobipen.wordpress.com/2018/11/27/visualize-pharmacophore-in-rdkit-rdkit/) is also great, even in 2D.
Unrelated question: do you have an idea on how to integrate this kind of widget in a notebook (https://codepen.io/diallobakary4/pen/MZzKZz?editors=0010).
https://stackoverflow.com/questions/54157893/integrate-javascript-widget-in-jupyter-notebook