Draw a cloud or lines over the polar area of a molecule in RDKit

137 views Asked by At

cordial greeting. I'm more of a chemist than a programmer. I am trying to use RDKit to draw in highlight or in some way (as a program called openeye does) the polar surface area of a molecule in two dimensions. I have the following code, I don't know if this can be improved in any way, there should be no polar area on the aromatic ring, but the program still gives me red dots, I understand that it is due to the partial charge calculated by the Gasteiger model

from rdkit import Chem
from rdkit.Chem import Draw
from rdkit.Chem import AllChem
from rdkit.Chem import Descriptors
#from rdkit.Chem import DrawMol

def highlight_psa_atoms(mol):
    # Gasteiger charges
    AllChem.ComputeGasteigerCharges(mol)
    
    # 
    psa_atoms = [atom.GetIdx() for atom in mol.GetAtoms() if atom.GetDoubleProp("_GasteigerCharge") < 0]
    
    # 
    highlight_style = {atom_id: (1, 0, 0) for atom_id in psa_atoms}
    
    return highlight_style

# 
smiles = "CC(=O)OC1=CC=CC=C1C(O)=O"

# 
mol = Chem.MolFromSmiles(smiles)

# 
highlight_style = highlight_psa_atoms(mol)

# 
img = Draw.MolToImage(mol, size=(300, 300), highlightAtoms=highlight_style, wedgeBonds=True, kekulize=True, wedgeLineWidth=2)
img
2

There are 2 answers

2
Tarquinius On

I know the output is way less beautiful than it is with openeye but maybe that's what you are looking for.

I'd recommend to directly visualize the contributions to the TPSA. Because the ring is not contributing to the TPSA it has no highlighting.

By default atom species P and S are (in terms of RDKit) not contributing to the TPSA. You can turn this on and off with the parameter includeSandP=True.

from rdkit import Chem
from rdkit.Chem import Draw
from rdkit.Chem.Draw import rdMolDraw2D
from rdkit.Chem import rdMolDescriptors
from IPython.display import Image

# Load a molecule
smiles = "CC(=O)OC1=CC=CC=C1C(O)=O"
mol = Chem.MolFromSmiles(smiles)

# Calculate the TPSA and identify atoms contributing to TPSA
# Include S and P, or don't - as you wish - here
tpsa_contribs = rdMolDescriptors._CalcTPSAContribs(mol, includeSandP=True)

# Find atoms that contribute to TPSA (usually N and O atoms)
highlight_atoms = [i for i, contrib in enumerate(tpsa_contribs) if contrib > 0]

# Create a drawer object for PNG output
drawer = rdMolDraw2D.MolDraw2DCairo(300, 300)

# Draw the molecule with highlighted atoms
drawer.DrawMolecule(mol, highlightAtoms=highlight_atoms)
drawer.FinishDrawing()

# Get the PNG data from the drawer
png_data = drawer.GetDrawingText()

# Display the PNG in the Jupyter notebook
Image(png_data)

Good find, GetSimilarityMapFromWeights allows you to put weights on certain atoms. Depending on the weights you can color it. Best in your case would be a diverging colormap. You can find all the colormaps here. Try the following:

import numpy as np
from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, rdMolDescriptors
from rdkit.Chem.Draw import SimilarityMaps

# Load a molecule
smiles = "CCNC(=O)NC1=NC2=CC=C(C=C2S1)C(=O)NCCS" 
mol = Chem.MolFromSmiles(smiles)

# Calculate the TPSA
# Include S and P, or don't - as you wish - here
tpsa = rdMolDescriptors._CalcTPSAContribs(mol, includeSandP=True)

# Get the image
fig = SimilarityMaps.GetSimilarityMapFromWeights(
    mol,
    size=(400, 400),
    weights=tpsa,
    colorMap='bwr',  # play around
    contourLines=10  # play around
)
fig.savefig('tpsa.png',bbox_inches='tight')

Give it a try and let me know.

0
Time Step On

Thank you very much @Tarquinius, I found codes that complement your answer:

import numpy as np
import matplotlib.pyplot as plt
from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, rdMolDescriptors
from rdkit.Chem.Draw import SimilarityMaps
# Load a molecule
smiles = "CC(=O)OC1=CC=CC=C1C(O)=O" 
mol = Chem.MolFromSmiles(smiles)
# Calculate the TPSA
tpsa = rdMolDescriptors._CalcTPSAContribs(mol)
# Get the image
fig = SimilarityMaps.GetSimilarityMapFromWeights(mol,size=(150, 150),weights=tpsa)
fig.savefig('tpsa.png',bbox_inches='tight')

The result is here

Thanks in advance to @Tarquinius