アロマティックなヘテロ環を漏れなくさがす

MORPHというアロマティックなヘテロ環を探索するツールの論文があります。

結局、総当りで探していってバリデータで確認すればいいっぽいので、単環で探すようなものを書いてみました。

戦略としては単純で、まず総当りで評価するSMILES Stringを生成します。原子のタイプ(atomchoices)とボンドのタイプ(bondchoices)を別々に生成してからコンビネーションをとれば、全ての組み合わせがつくれます(もれなくだぶりありなのでミッシーではないですが)。

アロマティックなリングかどうかを判定するのに、アロマティックなリングは小文字になるという性質を使います。要するにopenbabelでSMILES文字列を読んで、SMILES文字列を出力させてみて小文字だったら、芳香環が生成できたと判断します。find_aromatic_ringsという部分です。

from pybel import *
import re

def atomchoices(n):
    if n == 1:
        return ["C","N","O"]
    else:
        return [atom+atoms for atoms in atomchoices(n-1) for atom in ["C","N","O"]]
       
def bondchoices(n):
    if n == 1:
        return ["-","="]
    else:
        return [bond+bonds for bonds in bondchoices(n-1) for bond in ["-","="]]

def make_smiles(n):
    atoms = atomchoices(n)
    bonds = bondchoices(n)
    atoms = [a for a in atoms if "C" in a]
    bonds = [b for b in bonds if b.count("=") / float(n) <= 0.5]

    for a in atoms:
        for b in bonds:
            smiles = ""
            for i in range(n):
                smiles = smiles + a[i]
                if i == 0: smiles += "1"
                smiles = smiles + b[i]               
                if i == n-1: smiles += "1"
            if "=O" not in smiles and "O=" not in smiles and "O-O" not in smiles:
                mymol = str(readstring("smi",smiles))[:-2]
                yield mymol

def find_aromatic_rings(n):
    rearoma = re.compile('^[a-z1-9]+$')
    aromatic_compounds = list(set([smi for smi in make_smiles(n) if rearoma.search(smi)]))
    print aromatic_compounds

縮合環でやるならそれ用のテンプレート文字列を渡せるようにするか内部で辞書的に持っておいてそれを置き換えるようにすればいいんじゃないかと思います。