openbabelのPythonのとこ見てたらFingerprintsのとこで

>>> print fps[0] | fps[1] # Print the Tanimoto coefficient
0.3333

タニモト距離を求めるように、メソッドオーバーライドしてんのかと。コードを眺めた。

class Fingerprint(object):
    def __init__(self, fingerprint):
        self.fp = fingerprint
    def __or__(self, other):
        return ob.OBFingerprint.Tanimoto(self.fp, other.fp)
    @property 
    def bits(self):
            return _findbits(self.fp, ob.OBFingerprint.Getbitsperint())
    def __str__(self):
        return ", ".join([str(x) for x in self.fp])

参考

• 3.4.7 数値型をエミュレーションする

Python クックブック 第2版
Alex Martelli,Anna Martelli Ravenscroft,David Ascher
オライリー・ジャパン / ￥ 4,410 ()
在庫あり。

openbabelのringsearchがちょっと使いにくい気がする。部分構造をmolオブジェクトで返すメソッドがあっても良いんではないか？

from openbabel import *

obconversion = OBConversion()
obconversion.SetInFormat("smi")
obconversion.SetOutFormat("smi")
obmol = OBMol()

notatend = obconversion.ReadString(obmol,"c1ccccc1CC(=O)C2CCC2C(=O)O")

for r in obmol.GetSSSR():
    for a in  r._path:
        print a

これだと原子のIDしか返さない。

でも縮環チェックしたいんだったら、原子のリストがかぶるかどうかをリングの構成原子のIDをつかって調べればいいのでminimumの環からmaximumの環にするのは難しくはないか。

アトムのリストを渡したら部分構造検索を実行して、結果としてマッチした構造をmolとして返すメソッドがあればいいのか。

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の出力ファイルから構造を構築するものを用意すれば、頻出する部分構造を扱えるようになる。

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

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

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

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

XYZでreadした分子(CC)もAddHydrogensできない

>>> from openbabel import *
>>> conv = OBConversion()
>>> conv.SetInFormat("xyz")
True
>>> mol = OBMol()
>>> conv.ReadFile(mol, "test.xyz")
True
>>> mol.NumAtoms()
2  
>>> mol.AddHydrogens()
True
>>> mol.NumAtoms()
2

でも分子中のそれぞれの原子を指定して

mol.AddHydrogens(atom)

としてやればきちんと水素が付加されることがわかったので、とりあえずスクリプト中でループをまわせばいいかな

pythonラッパーで。openbabelで気をつけないといけないところはAtomのインデックスは1から始まるのに対しBondのインデックスは0からはじまる。

from openbabel import *

mol = OBMol()

a = mol.NewAtom()
a.SetAtomicNum(6)
a.SetVector(0.0, 1.0, 2.0)

b = mol.NewAtom()
b.SetAtomicNum(6)
b.SetVector(1.0, 2.0, 2.0)

mol.AddBond(1, 2, 2) 
mol.SetTitle("Ethene")

obConversion = OBConversion()
obConversion.SetOutFormat("smi")

smiles_string = obConversion.WriteString(mol)

print smiles_string

あとAddHydrogensで水素を付与した場合の水素原子はSetVectorメソッドで座標が貼付けられないの。よくわからん。サンプルだときちんと水素がついてるので、なんかの属性をセットしとく必要があるのではなかろうか。

Open Babel 2.2.0にバージョンがあがったのでインストールをした。ソースから。

同時に、perl,python,rubyのバインディングもコンパイルしてインストールしておく。

WindowsにおいてあるPlutoのファイル群を持ってくる

Windowsで開発しているPlutoをmacでもいじれるようにする。windowsのほうはhg serveと打てばhttpサーバーが立ち上がり、port8000番でアクセスできるようになる。

macのほうはディレクトリを作って初期化してpullしてupdate

mkdir Pluto
cd Pluto
hg init
hg pull http://192.168.XXX.XXX:8000/
hg update

これでOK