rdkit / rdkit

The official sources for the RDKit library
BSD 3-Clause "New" or "Revised" License
2.67k stars 881 forks source link

GetAcceptor1FeatVects in FeatDirUtilsRD reverses conditional cases (e.g. for carbonyls and fluorines) #3433

Open kangway opened 4 years ago

kangway commented 4 years ago

RDKit Version: '2020.03.3' Operating system: MacOS, Linux Are you using conda? yes If you are using conda, which channel did you install the rdkit from? conda-forge

I wanted to quickly test out the pharmacophore features as shown in this blog post and noticed that the encoded feature directions don't seem to point in the right direction.

I tested this out on lorlatinib, which has both cases for it's aromatic fluorine and carbonyl. E.g. from the Chem.Features.FeatDirUtils.py:

def GetAcceptor1FeatVects(conf, featAtoms, scale=1.5) :
  """
  Get the direction vectors for Acceptor of type 1
  This is a acceptor with one heavy atom neighbor. There are two possibilities we will
  consider here
  1. The bond to the heavy atom is a single bond e.g. CO
     In this case we don't know the exact direction and we just use the inversion of this bond direction
     and mark this direction as a 'cone'
  2. The bond to the heavy atom is a double bond e.g. C=O
     In this case the we have two possible direction except in some special cases e.g. SO2
     where again we will use bond direction

As shown on the right, the aromatic fluorine gets tagged as case 2 instead of case 1 and the carbonyl vice versa (i.e. we want the fluorine to have the directional cone, and the carbonyl to have two possible directions for their lone pairs):

image

I believe the reversal is caused by line 380 of FeatDirUtilsRD.py, where '>' should be replaced with '==' to satisfy the condition on line 385:

singleBnd = mol.GetBondBetweenAtoms(aid,heavyAt.GetIdx()).GetBondType() > Chem.BondType.SINGLE

  # special scale - if the heavy atom is a sulfur (we should proabably check phosphorous as well)
  sulfur = heavyAt.GetAtomicNum()==16

  if singleBnd or sulfur:
    v1 = conf.GetAtomPosition(heavyAt.GetIdx())

Updating this line accordingly seems to fix this issue:

singleBnd = mol.GetBondBetweenAtoms(aid,heavyAt.GetIdx()).GetBondType() == Chem.BondType.SINGLE

image

h4RkhW8t53e commented 3 years ago

Hi, Looking at your diagram, I was also wondering if:

  1. The NH2 should perhaps have two vectors in the directions of the H atoms?
  2. Is the amide N an acceptor?
github-actions[bot] commented 3 days ago

This issue was marked as stale because it has been open for 90 days with no activity.