pythonimagerdkit

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.

enter image description here


Solution

  • 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:

    enter image description here

    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:

    enter image description here

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