SDFからGASTON実行できるようにしたので、結果を解釈するためにsmilesに変換するコードを用意した。

import re
import openbabel as ob

def convert(gf):
    obc = ob.OBConversion()
    obc.SetOutFormat('smi')
    mol = ob.OBMol()

    for l in gf.split('\n'):
        if len(l) > 0 and l[0] == 'v':
            a = mol.NewAtom()
            atomic_num = int(l.split(' ')[2])
            a.SetAtomicNum(atomic_num)

        elif len(l) > 0 and l[0] == 'e':
            begin_atom_idx = int(l.split(' ')[1]) + 1
            end_atom_idx   = int(l.split(' ')[2]) + 1
            bond_order     = int(l.split(' ')[3])
            b = mol.AddBond(begin_atom_idx, end_atom_idx, bond_order)

        elif len(l) > 0 and l[0] == '#':
            title = l.split(' ')[1]
            mol.SetTitle(title)

    return obc.WriteString(mol)

if __name__ == '__main__':
    txt = open('gaston.out').read()
    p = re.compile('#.+?(?=(#|$))',re.S)
    m = p.finditer(txt)

    for ss in m:
        print convert(ss.group())[:-1]

実行結果

NC  11
NCC 11
NCCC    11
NCC(=C)C    11
NCCC=C  11
NCC(=C)C=C  11
NCCC=CC 11
NCC(=C)C=CC 11
NCC=C   11

頻度の高いフラグメントのSMILESと、その分子のIDが出力される

それにしてもpythonの正規表現はややこしいなぁと、みんなのpythonを読み返して復習

みんなのPython
柴田 淳
ソフトバンククリエイティブ / ￥ 2,940 ()

正規表現文字列 -> コンパイル -> 正規表現オブジェクト -> マッチ処理 -> マッチオブジェクト

gSpanやGASTONの入力はラベル付きのvertexとedgeなので、openbabelでsdfを読み込んで、GASTONのinputに変換するものを作ってみた。

import openbabel as ob

def convert(sdf):
    obc = ob.OBConversion()
    obc.SetInAndOutFormats('sdf','smi')

    mol = ob.OBMol()
    next = obc.ReadFile(mol,sdf)
    molnum = 0
    while next:
        # mol
        print "t # %d" % molnum
        # atom
        for i,atom in enumerate(ob.OBMolAtomIter(mol)): 
            print "v %d %d " % (i,atom.GetAtomicNum())

        # bond
        for i,bond in enumerate(ob.OBMolBondIter(mol)): 
            print "e %d %d %d" % (bond.GetBeginAtomIdx()-1,bond.GetEndAtomIdx()-1,bond.GetBondOrder())

        mol = ob.OBMol()
        next = obc.Read(mol)

        molnum += 1
    return True

if __name__ == "__main__":
    sdffile = 'pubchem_sample.sdf'

    convert(sdffile)

SDFはpubchemから適当に選んだがCID: 16757835は外しといた。

./gaston 11 pubchem_data pubchem_out

GASTONを実行した出力の一部

# 14
t 163
v 0 6
v 1 6
v 2 1
e 0 1 2
e 1 2 1
# 12
t 164
v 0 6
v 1 6
v 2 1
v 3 1
e 0 1 2
e 0 3 1
e 1 2 1

アロマティックなボンドは3みたいに別のラベルを与える必要がある気はするが、、、

後は、GASTONの出力ファイルから構造を構築するものを用意すれば、頻出する部分構造を扱えるようになる。

今度の発表の掴みを考えてみた。

htsはアッセイという関数を引数にとりコンパウンドリストを返す高階関数と考えることができる。

def hts(assay,library):
    criteria = 50 # %inh
    return [c for c in library if assay(c) > criteria]

うーん、マニアすぎるか

def chemist(comp):
    return modify_compound(c)

これまたあれだな。

おー、このタイムライン分かりやすい。

インフラとしてはディスクリプターをどう発生させとくか（毎度毎度モデル構築のたびに発生させるのはだるいし、基本的に一回作っとけば良いのでコンパウンドをどっかに登録するタイミングでバックグラウンドでやればいい仕事）とデータへのアクセス方法（いまはRESTっぽいのがいいんじゃないかと思ってる）をきちんとやっておけば、僕らはモデル構築の試行錯誤にそのリソースのほとんどを投入できていいよね。

RESTful Webサービス
Leonard Richardson,Sam Ruby
オライリー・ジャパン / ￥ 3,990 ()
在庫あり。

あと、こういう状況だと、モデル構築能力的には常に最先端のアプローチを試せるような状況にしとかないとまずいから、自然とデフォルトの環境がRになるような気がするし、立ち位置もデータマイニング寄りになると思うんだよね。

そして、立ち位置がマイニング寄りになっちゃうと、ケモインフォマティクスとかバイオインフォマティクスってのは、新しいモデル化の方法を試すための一つのターゲットに過ぎなくて、プロジェクトのモデル化とか創薬ビジネスモデル構築とかのほうに興味が移っちゃうのだよなぁ。

というエントリをケモインフォクックブックにのせておきました。

AMMOSというinduced fitを考慮したドッキングツールがopen sourceででるそうです。

AMMOS: Automated Molecular Mechanics Optimization tool for in silico Screening

ダウンロードできるようになったら試してみよう。DUDみたいなドッキングシミュレーションのベンチマーク用のデータベースもそろってきてるし、そろそろ今時のインターフェースで使いやすいドッキングシミュレーションツールのプロジェクトなど立ち上がりそうな気もしますが(もうなにか動いてるのかな)。

openbabelだったらobgenをつかえばいい。

pythonバインディングからやりたかったので、obgen.cppをよんだらOBBuilderクラスのインスタンスを生成すればいいらしいのだが、pythonからは現状OBBuilderのインスタンスをつくれないらしい。

MLを検索すると似たようなやりとりがあって、ラッパーを使えばいいらしく、pybelのソース見ろやーって書いてあったのでみた。

OPSってのは新たに導入された簡易プラグインシステムらしい。

で、ケモインフォクックブックにコードのせといた。

bioinformaticsの学習道路はそこそこ充実している。

でも、chemoinformaticsだってオープンソースで色々できるし、独学でも結構いけるよ！なんて思うんだけど情報が足りないなぁと思った。

というわけで、ちょっと小道でも。

ケモインフォクックブック

特にperl,pythonでちょっとした化学情報をいじるコードをクックブック形式で書いていけたらなぁなんて思ってる。

週に一回くらいは更新するつもりで。

MolBlasterは細切れにして、頻度を返さないといけないんだが、Separateメソッドがバグッてるらしいので、canonical smilesで出してpythonのほうで.でsplitして頻度をみてみた

import openbabel as ob
from random import sample

def randomsplit(mol,cutnum=5):
    cutlist =  sample(range(mol.NumBonds()),cutnum)

    delbonds = []
    for i,bond in enumerate(ob.OBMolBondIter(mol)): 
        if i in cutlist: delbonds.append(bond)
    for b in delbonds: mol.DeleteBond(b)

    return mol


def molblast(smi,iter=100,cutnum=5):
    obc = ob.OBConversion()
    obc.SetInAndOutFormats('smi','smi')
    obc.AddOption('c',ob.OBConversion.OUTOPTIONS)
    freq = {}

    for i in range(iter):
        mol = ob.OBMol()       
        obc.ReadString(mol,smi)
        for fragment in obc.WriteString(randomsplit(mol,cutnum=cutnum))[:-2].split('.'):
            freq[fragment] = freq.get(fragment,0) + 1

    return freq

if __name__ == "__main__":
    smiles = 'CCC(C)C(=O)OC1CC(C=C2C1C(C(C=C2)C)CCC(CC(CC(=O)O)O)O)O'

    print molblast(smiles,iter=1000,cutnum=20)

こんな感じ。

{'CCCCC=CC': 2, 'O=COCCC': 1, 'CCC=C(C)C(C)CC': 1, 'CCC(=O)O': 26, 'CCCCC(C)O':
 2, 'CCO': 459, 'CC=CCC(C)C': 2, 'CCC=CC(C)O': 2, 'CCCCC(C)C': 4, 'CCC': 705, 
'CCC=CC=CCC': 1, 'OCCCCC': 10, 'OCCC(C)OC': 1, 'C=CC=CCC': 2, 'CCCC(=C)C': 3, 
...

細切れにしないと頻度が稼げないし、細切れにするとファーマコフォア的な構造特徴が失われてしまう感じだ。というわけで、やってみて分かったけど、多分この方法だと実用的じゃないなぁ。Fragment Dependency Graphとかは面白そうなんだけどなぁ。

Molblasterっていうのは化学構造をランダムに切断していくっていう試行を何度か繰り返す方法で、その後フラグメントの頻度をシャノンエントロピーで評価したりするらしい。

シャノンエントロピーのほうに興味があったのとちょっと簡単に実装できそうだったのでpython+openbabelでやってみた。

構造をランダムに切断したい

まぁでも、構造活性相関にもっていくんだったら普通はベイズの方にいくよなぁ。多分著者の論文あると思うんだけど調べてない。