Dear James,
On Thu, Jun 10, 2010 at 2:35 PM, James Davidson <j.david...@vernalis.com> wrote:
>
> I have been trying figure-out how to return the count of aromatic rings for
> molecules (in Python), and am going to have to admit defeat! I saw in an
> earlier message
> (http://www.mail-archive.com/rdkit-discuss@lists.sourceforge.net/msg00153.html)
> a similar query, but I'm afraid it didn't help me very much. I also read
> the section on Aromaticity in the rdkit book, and realised that maybe this
> isn't a trivial exercise!
Correct. Counting the number of non-fused rings that are aromatic,
like the post you reference does, is pretty easy; including the fused
rings that are aromatic is more challenging.
> I would like the count to count aromatic ring-systems such that bicyclic (eg
> indole or naphthalene) would only count as 1. For reference, this appears
> to be the behaviour of the OpenEye OEDetermineAromaticRingSystems function -
> where the molecule derived from the smiles
> "C(O)(=O)c1cccc2c1[nH]c(C3CCCc4c3cccc4)c2" (which contains an indole and a
> tetrahydronaphthalene) gives a count of 2.
>
> Any help would be greatly appreciated.
I've attached a script that's not quite what you want, but it gets you
almost there: it finds all aromatic ring systems, including fused
ones. Anthracene, for example, gives 6 rings. The modifications to
this to get what you're looking for aren't a straightforward
post-processing step, but shouldn't be too bad. If there's not enough
here, let me know and I will take a look at adding the extra code.
This code isn't perfectly polished and could certainly be faster, but
it does seem mostly functional.
-greg
#
# Copyright (C) 2010 Greg Landrum
#
from rdkit import Chem
def IsRingAromatic(ring,aromaticBonds):
for bidx in ring:
if not aromaticBonds[bidx]:
return False
return True
def HasRingAromatic(ring,aromaticBonds):
for bidx in ring:
if aromaticBonds[bidx]:
return True
return False
def GetFusedRings(rings):
res=rings[:]
pool=[]
for i in range(len(rings)):
for j in range(i+1,len(rings)):
ovl=rings[i]&rings[j]
if ovl:
fused=rings[i]|rings[j]
fused.difference_update(ovl)
pool.append(fused)
while pool:
res.extend(pool)
nextRound=[]
for ringi in rings:
li=len(ringi)
for poolj in pool:
ovl = ringi&poolj
if ovl:
lj = len(poolj)
fused = ringi|poolj
fused.difference_update(ovl)
lf = len(fused)
if lf>li and lf>lj and fused not in nextRound and fused not in res:
nextRound.append(fused)
pool = nextRound
return res
def FindAromaticRings(mol):
# flag whether or not bonds are aromatic:
aromaticBonds = [0]*mol.GetNumBonds()
for bond in mol.GetBonds():
if bond.GetIsAromatic():
aromaticBonds[bond.GetIdx()]=1
# get the list of all rings:
ri = mol.GetRingInfo()
# collect the ones that have at least one aromatic bond:
rings=[set(x) for x in ri.BondRings() if HasRingAromatic(x,aromaticBonds)]
# generate all fused ring systems from that set
fusedRings=GetFusedRings(rings)
aromaticRings = [x for x in fusedRings if IsRingAromatic(x,aromaticBonds)]
return aromaticRings
if __name__=='__main__':
m = Chem.MolFromSmiles('C1=CC2=CC=CC=CC2=C1')
print FindAromaticRings(m)
m = Chem.MolFromSmiles('C1=CC2=C(C=C1)C=C1C=C3C=C4C=C5C(C=CC=C5C5=CC=CC6=C5C=CC=C6)=CC4=CC3=CC1=C2')
print FindAromaticRings(m)
------------------------------------------------------------------------------
ThinkGeek and WIRED's GeekDad team up for the Ultimate
GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the
lucky parental unit. See the prize list and enter to win:
http://p.sf.net/sfu/thinkgeek-promo
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
