Unexpected behavior of do_enhanced_stereo_sss?

[someonetookmyname]

Hi,

My problem is this: I enumerate the 2 isomers of a smiles having a single "and" group, and then I search against the results using one of the isomers.

My ability to find a match seems to require do_enahanced_stereo_sss = false, but I'm not sure why that would be, since the enumerations have no enhanced stereochemistry in them (also since I'm doing an exact match, not a substructure)?

The below example recreates the behavior .. apologies that it involves both API and cartridge .. that's the context I'm developing under. Which I've also retained in its Flask context, just so I can be sure the example code is valid. I'm using 2023_09_4.

The below code generates this response:

sudo curl --header "Content-Type: application/json" --request POST http://localhost:5001/testcase
{
  "enumerated": [
    "Cl[C@H]1CC[C@@H](Br)CC1",
    "Cl[C@H]1CC[C@H](Br)CC1"
  ],
  "pksfalse": [1],
  "pkstrue": []
}

Thank you for any information.

Best,

Andrew

@app.route('/testcase', methods = ['POST'])
def testcase():
    toret = {}
    pkstrue = []
    pksfalse = []
    enumerated = []

    conn = psycopg2.connect("abc123")
    curs = conn.cursor()
    curs.execute('drop table if exists tst')
    curs.execute('create table tst (pk integer, m mol)')
    curs.execute('create index tstidx on tst using gist(m)')

    # read-in/prepare the smiles
    _sp = Chem.SmilesParserParams()
    _sp.allowCXSMILES = True
    _sp.debugParse = False
    _sp.parseName = False
    _sp.removeHs = False
    _sp.sanitize = False
    _sp.strictCXSMILES = True
    m = Chem.MolFromSmiles("Cl[C@H]1CC[C@H](Br)CC1 |a:4,&1:1|", _sp)
    AllChem.Compute2DCoords(m)
    Chem.SanitizeMol(m)

    # enumerate the smiles molecule
    se_options = StereoEnumerationOptions(tryEmbedding=False, onlyUnassigned=True, maxIsomers=0, rand=None, unique=True, onlyStereoGroups=True)
    isomers = tuple(EnumerateStereoisomers(m, options=se_options))
    _sparams = Chem.SmilesWriteParams()
    _sparams.doRandom = False
    _sparams.rootedAtAtom = -1
    _sparams.doIsomericSmiles = True
    _sparams.canonical = False
    _sparams.doKekule = False
    _sparams.allHsexplicit = False
    _sparams.allBondsExplicit = False

    # store the enumerations
    for idx, mol in enumerate(isomers):
        curs.execute('insert into tst (pk, m) values (%s, mol_from_pkl(%s))', (idx, mol.ToBinary()))
        enumerated.append(Chem.MolToCXSmiles(mol, _sparams, flags=(CXSmilesFields.CX_ALL ^ CXSmilesFields.CX_COORDS)))
    conn.commit()

    curs.execute("set rdkit.do_chiral_sss=true")

    # search for one of the enumerations under do_enhanced_stereo_sss=true
    curs.execute("set rdkit.do_enhanced_stereo_sss=true")
    curs.execute('select pk from tst where m @= %s', (enumerated[1], ))
    res = curs.fetchall()
    for r in res:
        pkstrue.append(r[0])

    # search for one of the enumerations under do_enhanced_stereo_sss=false
    curs.execute("set rdkit.do_enhanced_stereo_sss=false")
    curs.execute('select pk from tst where m @= %s', (enumerated[1], ))
    res = curs.fetchall()
    for r in res:
        pksfalse.append(r[0])
    curs.close();
    conn.close();
    toret['pkstrue'] = pkstrue
    toret['pksfalse'] = pksfalse
    toret['enumerated'] = enumerated
    return toret

[greglandrum]

@someonetookmyname this is too much code for me to wade through in order to try and figure out what is going on.

If you'd like help, please make it possible for us to do so by providing a simple reproducible which uses either Python or the cartridge.

[someonetookmyname]

Hi,

The below example looks to reproduce the issue: it appears that if the cxsmiles used to query the db has atom coordinates, then the query will return nothing unless do_enhanced_stereo_sss==false.

It is unclear to me why this dependency would exist, especially since the db and query structure are identical, neither uses enhanced stereochemistry, and it is an equality match.

Thank you for the help and insights.

Andrew

drop table if exists tst;
DROP TABLE
create table tst (pk integer, m mol);
CREATE TABLE
create index tstidx on tst using gist(m);
CREATE INDEX
insert into tst (pk, m) values (1, mol_from_ctab('
     RDKit          2D

  8  8  0  0  0  0  0  0  0  0999 V2000
    3.0000    0.0000    0.0000 Cl  0  0  0  0  0  0  0  0  0  0  0  0
    1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.0000    0.0000    0.0000 Br  0  0  0  0  0  0  0  0  0  0  0  0
   -0.7500    1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7500    1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  2  1  1  1
  2  3  1  0
  3  4  1  0
  4  5  1  0
  5  6  1  6
  5  7  1  0
  7  8  1  0
  8  2  1  0
M  END'));
INSERT 0 1
set rdkit.do_chiral_sss=true;
SET
set rdkit.do_enhanced_stereo_sss=true;
SET
select pk from tst where m @= 'Cl[C@H]1CC[C@H](Br)CC1'::mol;
 pk 
----
  1
(1 row)

select pk from tst where m @= 'Cl[C@H]1CC[C@H](Br)CC1 |(3,0,;1.5,0,;0.75,-1.29904,;-0.75,-1.29904,;-1.5,0,;-3,0,;-0.75,1.29904,;0.75,1.29904,)|'::mol;
 pk 
----
(0 rows)

set rdkit.do_enhanced_stereo_sss=false;
SET
select pk from tst where m @= 'Cl[C@H]1CC[C@H](Br)CC1'::mol;
 pk 
----
  1
(1 row)

select pk from tst where m @= 'Cl[C@H]1CC[C@H](Br)CC1 |(3,0,;1.5,0,;0.75,-1.29904,;-0.75,-1.29904,;-1.5,0,;-3,0,;-0.75,1.29904,;0.75,1.29904,)|'::mol;
 pk 
----
  1
(1 row)

[someonetookmyname]

As some more due diligence on this issue, I noticed the following sentence on https://www.rdkit.org/docs/RDKit_Book.html for CXSMILES:

If 2D coordinates are present in the CXSMILES, atomic stereo can be set using wUorwD to create wedged bonds.

As the smiles I used in the db example of the last post were generated by MolToCXSmiles(), I needed to first call Chem.WedgeMolBonds() in order for the wU/wD markups to be generated in the CXSMILES.

However, the revised CXSMILES of
Cl[C@H]1CC[C@H](Br)CC1 |(3,0,;1.5,0,;0.75,-1.29904,;-0.75,-1.29904,;-1.5,0,;-3,0,;-0.75,1.29904,;0.75,1.29904,),wU:1.0wD:4.4|
behaved identically to the CXSMILES with coordinates but no wU/wD markups (required do_enhanced_stereo==false to find the result).

So I'm afraid I can still not understand the result/behavior.

Related to this, is there a preferred order for calling WedgeMolBonds() and AssignStereochemistry()/FindPotentialStereo()?

Thank you.

Best, Andrew

[greglandrum]

I'd characterize this as a bug. Here's what's happening: the molecule equality code goes through a bunch of increasingly expensive checks to detect that the two molecules are different. The last of those steps is to generate a SMILES or CXSMILES (in the case of enhanced stereo) for each and compare them. Since CXSMILES includes coordinates by default, this is going to result in molecules with different atomic positions comparing unequal.

In general: @= is there because we had to provide a molecular equality operator, but I think it's generally better to do molecular comparison with SMILES (or CXSMILES if you want enhanced stereo) if you can.

[greglandrum]

@someonetookmyname it would be cool if you could create an issue for this so that you get credit for it.

[someonetookmyname]

Thanks for this, Greg. Yes, I'll create an issue, though I certainly don't need any credit .. smile

Based on your reply, what is then the preferred way to do a stereochemistry-less exact match involving molecules with CXSMILES-requiring features (i.e. the only difference between molecules is stereochemistry)? It would seem pretty straightforward to manually strip out @'s and double-bond-orientation indicators from a normal SMILES, but CXSMILES can also have wedging info etc., which makes me more hesitant towards manual manipulation to culminate in a match of CXSMILES.

Or is resorting to bi-directional stereochemistry-less substructure matches the way to go?

Thanks for any thoughts.

[someonetookmyname]

And just so it doesn't appear like I'm ignoring the obvious, perhaps @= is preferred when do_chiral_sss==false? Though, do_chiral_sss still appears to require double-bond geometry to be exact?

[someonetookmyname]

With some more testing, do_chiral_sss does appear to impact inclusion/exclusion of considering double-bond geometry in matches