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

2 thoughts on “Visualize pharmacophore with RDKit #RDKit #Pymol #Chemoinformatics

  1. 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.

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 )

Google photo

You are commenting using your Google 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.