I am attempting to reduce a molecule (attached SDF) to just its ring
systems using Code Block 1 at the bottom.
The problem is that when I get through the loops removing non-ring
atoms/bonds, and convert the EditableMol back to a Mol, I end up with 7
disjoint sets of atoms:
((0, 1, 2, 3, 4, 5, 6, 7, 8), (9, 10, 11, 12, 13, 14), (15,), (16,), (17,
20), (18,), (19,))
It appears that when I remove an atom from the EditableMol by index, the
indices are reassigned. I tried to test this with the inelegant code in
Code Block 2, which gives me the expected sets of atom indices with respect
to number and size:
((0, 1, 2, 3, 4, 5, 6, 7, 8), (9, 10, 11, 12, 13, 14), (15, 16, 17, 18, 19,
20))
- but it still fails to sanitize when I convert back to a Mol.
What am I missing here? Also, is there an easier (ie, existing) way to do
this? I'm just looking to reduce the molecule to its ring systems and
write those to an SD file.
-Kirk
Code block 1:
from rdkit import Chem
sdin = Chem.SDMolSupplier('test.sdf')
sdout = Chem.SDWriter('rings.sdf')
for m in sdin:
print len(m.GetBonds()),len(m.GetAtoms())
em = Chem.EditableMol(m)
for a in m.GetAtoms():
if ( not a.IsInRing() ):
em.RemoveAtom(a.GetIdx())
print a.GetIdx(), m.GetAtomWithIdx(a.GetIdx()).GetSymbol()
for b in m.GetBonds():
if ( not b.IsInRing() ):
a1 = b.GetBeginAtomIdx()
a2 = b.GetEndAtomIdx()
em.RemoveBond(a1,a2)
m3 = em.GetMol()
print len(m3.GetBonds()), len(m3.GetAtoms())
f = Chem.GetMolFrags(em.GetMol())
print f
#for f in Chem.GetMolFrags(m3,asMols = True):
#sdout.write(f)
Code block 2:
from rdkit import Chem
sdin = Chem.SDMolSupplier('test.sdf')
sdout = Chem.SDWriter('rings.sdf')
for m in sdin:
print len(m.GetBonds()),len(m.GetAtoms())
em = Chem.EditableMol(m)
active = True
while ( active == True ):
active = False
for a in m.GetAtoms():
if ( not a.IsInRing() ):
print a.GetIdx(), m.GetAtomWithIdx(a.GetIdx()).GetSymbol()
em.RemoveAtom(a.GetIdx())
active = True
m=em.GetMol()
em = Chem.EditableMol(m)
break
active = True
while ( active == True):
active = False
for b in m.GetBonds():
if ( not b.IsInRing() ):
active = True
a1 = b.GetBeginAtomIdx()
a2 = b.GetEndAtomIdx()
em.RemoveBond(a1,a2)
m=em.GetMol()
em = Chem.EditableMol(m)
break
m3 = em.GetMol()
print len(m3.GetBonds()), len(m3.GetAtoms())
f = Chem.GetMolFrags(em.GetMol())
print f
#for f in Chem.GetMolFrags(m3,asMols = True):
#sdout.write(f)
------------------------------------------------------------------------------
Get 100% visibility into Java/.NET code with AppDynamics Lite
It's a free troubleshooting tool designed for production
Get down to code-level detail for bottlenecks, with <2% overhead.
Download for free and get started troubleshooting in minutes.
http://p.sf.net/sfu/appdyn_d2d_ap2
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
