Maybe the logic is that sophisticated after all - Code Block 4 below.
This should satisfy the majority of cases, but probably not all. I'm open
to better approaches.
Code Block 4:
from rdkit import Chem
from rdkit.Chem import AllChem
sdin = Chem.SDMolSupplier('test.sdf')
sdout = Chem.SDWriter('rings.sdf')
for m in sdin:
em = Chem.EditableMol(Chem.Mol())
indexmap = {}
for a in m.GetAtoms():
if ( a.IsInRing() ):
indexmap[a.GetIdx()] = em.AddAtom(Chem.Atom(a.GetAtomicNum()))
for b in m.GetBonds():
if ( b.IsInRing() ):
em.AddBond(
indexmap[b.GetBeginAtomIdx()],indexmap[b.GetEndAtomIdx()],b.GetBondType() )
for a in m.GetAtoms():
if ( a.IsInRing() and a.GetAtomicNum() == 7 ):
for b in a.GetBonds():
if ( not b.IsInRing() ):
em.AddBond(em.AddAtom(Chem.Atom(1)), indexmap[a.GetIdx()],
Chem.BondType.SINGLE)
for nm in Chem.GetMolFrags(em.GetMol(), asMols=True):
AllChem.Compute2DCoords(nm)
sdout.write(nm)
On Fri, May 31, 2013 at 2:41 PM, Robert DeLisle <[email protected]> wrote:
> 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
