How to obtain fingerprint bits as SMARTS pattern in RDKit

Question

I am using RDKit to generate Morgan Fingerprints (similar to ECFP) and then obtaining the bit information. I need the bit information in order to generate a statistics of substructures for each position of fingerprint.

I am using the code that is given in the RDkit getting started page (here). My code looks something like this (unimportant parts are edited out):

from rdkit import Chem
from rdkit.Chem.rdMolDescriptors import GetHashedMorganFingerprint

# for loop begins
mol = Chem.MolFromSmiles('CCCNCCC=O') # or some other SMILES
bit_info = {}
fp = GetHashedMorganFingerprint(mol, 3, bitInfo=bit_info) # eqv to ECFP6

atomidx, radius = bit_info[118][0] # 118 is the position that I am studying for example
env = Chem.FindAtomEnvironmentOfRadiusN(mol, radius, atomidx) # returns rdkit vector
submol = Chem.PathToSubmol(mol,env)
# store SMARTS somehow

Now this submol should contain the substructure fragment that is hashed into the 118 bit position in fingerprint.

I want to convert this substructure fragment into a SMARTS pattern (not a SMILES string). RDkit does allow me to use Chem.MolToSmarts() on submol. However, if I run it on position 118, then I get this SMARTS, visualised by SMARTS plus:

But I want the SMARTS to include information on the bonds that are attached to the outside of the fragment, because that information is also used by the fingerprint. So I want my SMARTS to show this:

I know that these images can be produced by rdkit.Chem.Draw.DrawMorganBit:

But I cannot use images as I need to sort through ~ 3000 molecules. So I need those substructures in SMARTS form so that I can compare them.

Another big problem is that PathToSubmol() does not seem to work if you set radius=0 which is true for many bit positions in fingerprint. But DrawMorganBit can handle the cases where radius=0.

So, I essentially need what the DrawMorganBit is showing, but in SMARTS format.

After I have the SMARTS, I can essentially compare fragments and different bit positions (that's a part of analysis for my research). I cannot use only the submol because it has missing connectivity information which is relevant during the hashing, so it determines the bit position where the fragment appears.

How can I achieve this?

Answer 1

I'm not entirely sure I understand your question. You want a SMARTS that matches DrawMorganBit but with an extra * added?

If you look at the code for DrawMorganBit you can see it boils down to an internal method _getMorganEnv:

  if not mol.GetNumConformers():
    rdDepictor.Compute2DCoords(mol)
  bitPath = Chem.FindAtomEnvironmentOfRadiusN(mol, radius, atomId)

  # get the atoms for highlighting
  atomsToUse = set((atomId, ))
  for b in bitPath:
    atomsToUse.add(mol.GetBondWithIdx(b).GetBeginAtomIdx())
    atomsToUse.add(mol.GetBondWithIdx(b).GetEndAtomIdx())

  #  enlarge the environment by one further bond
  enlargedEnv = set()
  for atom in atomsToUse:
    a = mol.GetAtomWithIdx(atom)
    for b in a.GetBonds():
      bidx = b.GetIdx()
      if bidx not in bitPath:
        enlargedEnv.add(bidx)
  enlargedEnv = list(enlargedEnv)
  enlargedEnv += bitPath
  # … (more code to get coordinates, etc.)

In other words, you can do exactly the same thing to enlarge your submol if that's what you want.

Answer 2

Based on your code and the code found by Geoff.

Additionally you have to find the neighbors and change their identity.

from rdkit import Chem
from rdkit.Chem.rdMolDescriptors import GetHashedMorganFingerprint

mol = Chem.MolFromSmiles('CCCNCCC=O')

bit_info = {}
fp = GetHashedMorganFingerprint(mol, 3, bitInfo=bit_info)

atomidx, radius = bit_info[118][0]
env = Chem.FindAtomEnvironmentOfRadiusN(mol, radius, atomidx)

atomsToUse = set((atomidx, ))
for b in env:
    atomsToUse.add(mol.GetBondWithIdx(b).GetBeginAtomIdx())
    atomsToUse.add(mol.GetBondWithIdx(b).GetEndAtomIdx())

enlargedEnv = set()
for atom in atomsToUse:
    a = mol.GetAtomWithIdx(atom)
    for b in a.GetBonds():
      bidx = b.GetIdx()
      if bidx not in env:
        enlargedEnv.add(bidx)

enlargedEnv = list(enlargedEnv)
enlargedEnv += env

# find all relevant neighbors
anyAtoms = []
for a in atomsToUse:
    neighbors = mol.GetAtomWithIdx(a).GetNeighbors()
    for n in neighbors:
        anyIdx = n.GetIdx()
        if anyIdx not in atomsToUse:
            anyAtoms.append(anyIdx)

# replace atomic number to zero (there is no number for any atom)
for aA in anyAtoms:
    mol.GetAtomWithIdx(aA).SetAtomicNum(0)

submol = Chem.PathToSubmol(mol, enlargedEnv)

# change [0] to *
MorganBitSmarts = Chem.MolToSmarts(submol).replace('[#0]', '*')
print(MorganBitSmarts)

This gives you *-[#6]-[#6]-[#6]=[#8]