I just want to share my script, which I use for enumeration of
stereoisomers. Enumeration of double bonds was added quite recently and
thus I didn't test it extensively.
I put it on github: https://github.com/DrrDom/rdk
It seems to work well on quite complex queries like
CC1CCC(CC1)C1CC=C(C\C(=C/I)C(=CF)C1C1C[C@@H](C)CC=C1)C1CC[C@H](O)C(Br)C1
to generate all possible stereoisomers:
gen_stereo_rdkit3.py -i input.smi -o output.smi -t -d
to discard compounds with more than 4 undefined centers/double bonds:
gen_stereo_rdkit3.py -i input.smi -o output.smi -t -d -u 4
Try it on your own risk :)
If you will find mistakes let me know.
Pavel.
On 12/10/2016 03:44 AM, Brian Cole wrote:
So I may need a little guidance on this one from someone with a little
more historical knowledge of RDKit.
I found the findPotentialStereoBonds function inside Chirality.cpp
that appears to do what we want, detect double bonds with unspecified
E/Z stereo chemistry. That's at least what the comment in the file
leads me to believe:
// Find bonds than can be cis/trans in a molecule and mark them as "any"
...
void findPotentialStereoBonds(ROMol &mol, bool cleanIt) {
But when you look at the code, it's not setting bonds to
Bond::STEREOANY, it's setting them to Bond::STEREONONE here:
...
// proceed only if we either want to clean the stereocode on
this bond
// or if none is set on it yet
if (cleanIt || dblBond->getStereo() == Bond::STEREONONE) {
dblBond->setStereo(Bond::STEREONONE);
...
Seems odd to me check that the stereo is Bond::STEREONONE, and then
set it to that value.
While findPotentialStereoBonds isn't exposed in Python, there is some
Java docs written for it in the Java wrappers:
%javamethodmodifiers RDKit::MolOps::findPotentialStereoBonds (
ROMol & mol, bool cleanIt = false )
"
/**
<p>
finds bonds that could be cis/trans in a molecule and mark them as
Bond::STEREONONE
<p>
So that documentation is actually correct. It does find bonds that
could be cis/trans and mark them as Bond::STEREONONE. That just isn't
very useful since all bonds are marked Bond::STEREONONE by default.
The only direct test case of findPotentialStereoBonds I can find is
the following:
MolOps::findPotentialStereoBonds(*m, false);
TEST_ASSERT(m->getNumAtoms() == 15);
TEST_ASSERT(m->getBondWithIdx(1)->getBondType() == Bond::DOUBLE);
TEST_ASSERT(m->getBondWithIdx(1)->getStereoAtoms().size() == 2);
TEST_ASSERT(m->getBondWithIdx(3)->getBondType() == Bond::DOUBLE);
TEST_ASSERT(m->getBondWithIdx(3)->getStereoAtoms().size() == 2);
That doesn't even test the getStereo() return value. Though it does
test getStereoAtoms() return value. So it looks like the proper thing
to do to detect unspecified bond stereo is to check for a non-empty
getStereoAtoms() vector?
So ideally, I would like to write an bond stereo enumerator that tests
if a double bond has unspecified stereo (testing the size of
getStereoAtoms() would work) and then sets the stereo, that is:
bond.SetStereo(Bond::STEREOE);
or
bond.SetStereo(Bond::STEREOZ);
And then hopefully the Chem.MolToSmiles will "do the right thing".
However, this commented out code in the wrapper has given me pause:
// this is no longer exposed because it requires that stereo atoms
// be set. This is a task that is tricky and "dangerous".
//.def("SetStereo",&Bond::setStereo,
// "Set the CIP-classification of the bond as a
BondStereo\n")
So apparently I shouldn't expose an "easy" way to do this. What is the
trickiness and dangerousness of this API? And could we make an easy
way to enumerate bond stereo?
Thanks!
On Fri, Dec 9, 2016 at 5:44 PM, Brian Cole <col...@gmail.com
<mailto:col...@gmail.com>> wrote:
This has me quite curious now, how do we detect unspecified bond
stereo chemistry in RDKit?
m = Chem.MolFromSmiles("FC=CF")
assert m.HasProp("_StereochemDone")
for bond in m.GetBonds():
print(bond.GetBondDir(), bond.GetStereo())
Yields:
(rdkit.Chem.rdchem.BondDir.NONE, rdkit.Chem.rdchem.BondStereo.STEREONONE)
(rdkit.Chem.rdchem.BondDir.NONE, rdkit.Chem.rdchem.BondStereo.STEREONONE)
(rdkit.Chem.rdchem.BondDir.NONE, rdkit.Chem.rdchem.BondStereo.STEREONONE)
Trying to force the issue:
Chem.AssignStereochemistry(m, cleanIt=True, force=True,
flagPossibleStereoCenters=True)
assert m.HasProp("_StereochemDone") for bond in m.GetBonds():
print(bond.GetBondDir(), bond.GetStereo())
Still yields:
(rdkit.Chem.rdchem.BondDir.NONE, rdkit.Chem.rdchem.BondStereo.STEREONONE)
(rdkit.Chem.rdchem.BondDir.NONE, rdkit.Chem.rdchem.BondStereo.STEREONONE)
(rdkit.Chem.rdchem.BondDir.NONE, rdkit.Chem.rdchem.BondStereo.STEREONONE)
I can see the setStereo(Bond::STEREOANY) in the code here:
if ((*bondIt)->getBondDir() == Bond::EITHERDOUBLE) {
(*bondIt)->setStereo(Bond::STEREOANY);
Guess I have to figure out how to detect EITHERDOUBLE?
-Brian
On Fri, Dec 9, 2016 at 5:03 PM, Andrew Dalke
<da...@dalkescientific.com <mailto:da...@dalkescientific.com>> wrote:
On Dec 9, 2016, at 9:50 PM, Brian Kelley wrote: > >>> from
rdkit import Chem > >>> m = Chem.MolFromSmiles("F/C=C/F") >
>>> for bond in m.GetBonds(): > ... print bond.GetStereo()
> ... > STEREONONE > STEREOE > STEREONONE > > However, setting
bond stereo doesn't appear to be exposed. I thought I knew a
way to change it, via GetBondDir()/SetBondDir() on the
directional bonds. It doesn't seem to work. First, get the
existing bond directions: >>> from rdkit import Chem >>> mol =
Chem.MolFromSmiles("F/C=C/F") >>> for bond in mol.GetBonds():
... print(bond.GetBondDir()) ... ENDUPRIGHT NONE ENDUPRIGHT
I'll change the first "/" to a "\" (That's "\\" when escaped
for a normal Python string.) >>> mol.GetBondWithIdx(0)
<rdkit.Chem.rdchem.Bond object at 0x10fe26600> >>>
mol.GetBondWithIdx(0).GetBondDir()
rdkit.Chem.rdchem.BondDir.ENDUPRIGHT >>>
mol.GetBondWithIdx(0).SetBondDir(Chem.BondDir.ENDDOWNRIGHT) If
I generate the SMILES I see the old "F/C=C/F" instead of the
new direction. I need to clear internal computed properties
when I make structure edits: >>> Chem.MolToSmiles(mol,
isomericSmiles=True) 'F/C=C/F' >>> mol.ClearComputedProps()
>>> Chem.MolToSmiles(mol, isomericSmiles=True) 'F/C=C\\F' I'll
set the directions to "NONE" and clear computed properties.
The SMILES is correct: >>>
mol.GetBondWithIdx(0).SetBondDir(Chem.BondDir.NONE) >>>
mol.GetBondWithIdx(2).SetBondDir(Chem.BondDir.NONE) >>>
mol.ClearComputedProps() >>> Chem.MolToSmiles(mol,
isomericSmiles=True) 'FC=CF' although the stereo setting is
unchanged: >>> for bond in mol.GetBonds(): ...
print(bond.GetStereo()) ... STEREONONE STEREOE STEREONONE I
expected it to give me all STEREONONEs, as if it were the same
as the following: >>> mol2 = Chem.MolFromSmiles("FC=CF") >>>
for bond in mol2.GetBonds(): ... print(bond.GetStereo()) ...
STEREONONE STEREONONE STEREONONE
Andrew da...@dalkescientific.com
<mailto:da...@dalkescientific.com>
------------------------------------------------------------------------------
Developer Access Program for Intel Xeon Phi Processors Access
to Intel Xeon Phi processor-based developer platforms. With
one year of Intel Parallel Studio XE. Training and support
from Colfax. Order your platform today.http://sdm.link/xeonphi
_______________________________________________ Rdkit-discuss
mailing list Rdkit-discuss@lists.sourceforge.net
<mailto:Rdkit-discuss@lists.sourceforge.net>
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
<https://lists.sourceforge.net/lists/listinfo/rdkit-discuss>
------------------------------------------------------------------------------
Developer Access Program for Intel Xeon Phi Processors
Access to Intel Xeon Phi processor-based developer platforms.
With one year of Intel Parallel Studio XE.
Training and support from Colfax.
Order your platform today.http://sdm.link/xeonphi
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
------------------------------------------------------------------------------
Developer Access Program for Intel Xeon Phi Processors
Access to Intel Xeon Phi processor-based developer platforms.
With one year of Intel Parallel Studio XE.
Training and support from Colfax.
Order your platform today.http://sdm.link/xeonphi
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
