Highlighting smiles features with rdkit

from rdkit import Chem
from rdkit.Chem import ChemicalFeatures
from rdkit import rdBase
from rdkit.RDPaths import RDDataDir
# from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit.Chem import AllChem
import os

print(rdBase.rdkitVersion)

    
fdefFile = os.path.join(RDDataDir,'BaseFeatures.fdef')
featFact = ChemicalFeatures.BuildFeatureFactory(fdefFile)



mols = [Chem.MolFromSmiles('C1=NC(=C2C(=N1)N(C=N2)C3C(C(C(O3)CO)O)O)N')]

featslists = [featFact.GetFeaturesForMol(mol) for mol in mols]



    
def drawp4core(mol, feats):
    atoms_list = {}
    for feat in feats:
        atom_ids = feat.GetAtomIds()
        feat_type = feat.GetType()
        atoms_list[feat_type] = atom_ids
    return Draw.MolsToGridImage([mol]*len(atoms_list), legends=list(atoms_list.keys()), highlightAtomLists=list(atoms_list.values()))

    
im = drawp4core(mols[0], featslists[0])

print(type(im))

im.show()

I am trying to highlight pharmacophoric features in small molecules using RDKIT. I have the following code (copied from ref, check it for the full reproducible code).

This works but all the highlights are in separate images like the example. I want them all superimposed in one hopefully in different colours.

I tried to figure out a way to highlight all the features as per your question. My first attempts resulted in highlighting the entire molecule.

I decided to go for a different color slices option, I didn't find a way to get the color linked to the legend though; result here and code below:

Code:

import os
import rdkit
from io import BytesIO
from PIL import Image
from rdkit.RDPaths import RDDataDir
from rdkit import Chem
from rdkit.Chem import AllChem 
from rdkit.Chem import Draw
from rdkit import Geometry
from cairosvg import svg2png

print('\n-------------------------------')
print('\n rdkit Version : ', rdkit.__version__)
print('\n-------------------------------')




def drawp4core(mol, feats):
    
    AllChem.Compute2DCoords(mol) ## needed to use atom_point = mol.GetConformer().GetAtomPosition(atom) -- # AllChem.Compute2DCoords(mol)
    
    # drawer = Draw.rdMolDraw2D.MolDraw2DCairo(800,800, -1,-1) # need SVG image
    
    drawer = Draw.rdMolDraw2D.MolDraw2DSVG(800,800, -1,-1) #
    

    print('len features : ', len(feats))
    
    for i in feats:
        
        # print(type(i) , i.GetFactory() , i.GetType() , i.GetFamily())
        
        print(i.GetType() ,'---', i.GetFamily())
              
    
    
    atoms_list = {}
    
    for feat in feats:
        
        atom_ids = feat.GetAtomIds()
        

        feat_type = feat.GetType()
        
        atoms_list[feat_type] = atom_ids
    
    
    drawer.DrawMolecule(mol, legend=str([i for i in (atoms_list.keys())]).strip("[]").replace("'", ""),
                               
                                highlightAtoms = [], highlightAtomColors=[])  
    
    for feat in feats:
        
        atom_ids = feat.GetAtomIds()
        
        for atom in atom_ids: 
            
            atom_point = mol.GetConformer().GetAtomPosition(atom)
            
            
            # print('\natom :' , atom , atom_point)
        
            if feat.GetType() == 'SingleAtomDonor' :
                
                print(Geometry.Point2D(atom_point.x, atom_point.y))
                
                drawer.SetColour((.7,0,.7, 0.5))
                
                Draw.MolDraw2D.DrawArc( drawer , Geometry.Point2D(atom_point.x, atom_point.y), 0.7, 0.0  , 60.00 ,
                                       
                                        False)
                
            if feat.GetType() == 'SingleAtomAcceptor' :
                 
                print(Geometry.Point2D(atom_point.x, atom_point.y))
                 
                drawer.SetColour((.7,.7,.7, 0.5))
                 
                Draw.MolDraw2D.DrawArc( drawer , Geometry.Point2D(atom_point.x, atom_point.y), 0.7 , 60.0  , 120.00 ,
                                        
                                        False)
                
                
            if feat.GetType() == 'Imidazole' :
                 
                print(Geometry.Point2D(atom_point.x, atom_point.y))
                 
                drawer.SetColour((1, 1 ,0 , 0.5))
                 
                Draw.MolDraw2D.DrawArc( drawer , Geometry.Point2D(atom_point.x, atom_point.y), 0.7, 120.0  , 180.00 ,
                                        
                                        False)
                
            if feat.GetType() == 'Arom5' :
                 
                print(Geometry.Point2D(atom_point.x, atom_point.y))
                 
                drawer.SetColour((0, 1 ,0 , 0.5))
                 
                Draw.MolDraw2D.DrawArc( drawer , Geometry.Point2D(atom_point.x, atom_point.y), 0.7, 180.0  , 240.00 ,
                                        
                                        False)
                
                
            if feat.GetType() == 'Arom6' :
                 
                print(Geometry.Point2D(atom_point.x, atom_point.y))
                 
                drawer.SetColour((1, .6 , .2 , 0.5))
                 
                Draw.MolDraw2D.DrawArc( drawer , Geometry.Point2D(atom_point.x, atom_point.y), 0.7, 240.0  , 300.00 ,
                                        
                                        False)
                
    drawer.DrawMolecule(mol, legend=str([i for i in (atoms_list.keys())]).strip("[]").replace("'", ""),
                                
                                 highlightAtoms = [], highlightAtomColors=[]) 
    
    
    
    drawer.FinishDrawing()

    svg_data = drawer.GetDrawingText()
    
    # with open('output.svg', 'w') as f:
            
    #     f.write(svg_data)
        
    png = svg2png(bytestring=svg_data)

    pil_img = Image.open(BytesIO(png)).convert('RGBA')
    # pil_img.save('output/pil.png')
    
    print(type(pil_img))  

    pil_img.show()
    




molecule = [Chem.MolFromSmiles('C1=NC(=C2C(=N1)N(C=N2)C3C(C(C(O3)CO)O)O)N')]

print(molecule , type(molecule))

print(molecule[0] , type(molecule[0]))


fdefFile = os.path.join(RDDataDir,'BaseFeatures.fdef')
featFact = Chem.ChemicalFeatures.BuildFeatureFactory(fdefFile)


featslists = [featFact.GetFeaturesForMol(mol) for mol in molecule]

for mol in molecule:

    mol_img = drawp4core(molecule[0], featslists[0])

First

drawer.DrawMolecule(mol, legend=str([i for i in (atoms_list.keys())]).strip("[]").replace("'", ""),
                                
                                 highlightAtoms = [], highlightAtomColors=[]) 

is needed because I needed a drawing to add the DrawArc pics to the image; the second one is to re-draw the molecule on top of the circles slices.

I wasn't able to find any way to add concentric circles like:

instead of slices. As I mentioned before, I didn't find a way to get the color linked to the legend though.