11 01 2013 chemoinformatics Python openbabel Tweet
高速文字列解析の流れで、タイミングよく統数研チャンネルの資料が公開されていて、最後のほうのスライド(65枚目)にグラフの類似度検索の実装としてgwtが紹介されていたので遊んでみた。
データベースは、ChEMBLから取ってきた。なぜpubchemじゃないかというと面倒なので略。
入力のフォーマットがgSpanとかGASTONと一緒だったので、昔書いたのを流用しようとしたんだが、昔書いたpythonスクリプトはどうもOBMolAtomIterあたりがバグっていて動かないようだ(OBMolBondIterも動かなかった。)
File "sdf2gsp", line 17, in convert for i,atom in enumerate(ob.OBMolAtomIter(mol)): File "/usr/local/Cellar/python/2.7.3/lib/python2.7/site-packages/openbabel.py", line 4499, in __init__ if not self.iter.__bool__(): AttributeError: '_OBMolAtomIter' object has no attribute '__bool__'
しょうがないのでiterateしない版で書き換えた。
#!/usr/bin/env python import sys import pybel import openbabel as ob def convert(sdf): mols = pybel.readfile('sdf',sdf) for n, mol in enumerate(mols, 1): print "t # %d" % n for i in range(mol.OBMol.NumAtoms()): print "v %d %d " % (i+1, mol.OBMol.GetAtom(i+1).GetAtomicNum()) for j in range(mol.OBMol.NumBonds()): b = mol.OBMol.GetBond(j) print "e %d %d %d" % (b.GetBeginAtomIdx(), b.GetEndAtomIdx(), b.GetBondOrder()) if __name__ == "__main__": sdffile = sys.argv[1] convert(sdffile)
ちなみにこの形式に名前ついているんだろうか？babelとかでサポートされないと使い勝手が悪いんだけどなぁ。
これで、ChEMBLの137万件をgspにコンバートしてみたけど遅すぎるので後でHaskellで書きなおす。sdf単に舐めるだけだからOBMolに変換する必要もないし。
できたchembl_14.gspからindexファイルをつくる。600M(137万化合物)で10分ぐらいかかった@MBA。
./gwt-build -iteration 2 chembl_14.gsp chembl14
検索はみんな大好きImatinibをgspにしたものを使った。
$ ./gwt-search -kthreshold 0.8 chembl14 query.gsp reading index file:chembl14 end reading start building rank dictionary end building rank dictionary start transferring labels end transferring labels id:0 1080597:0.809449 712949:0.838758 413038:0.847826 \ 394513:0.826087 391369:0.819028 179393:0.815124 \ 175534:0.899647 48523:0.804348 47164:1.000000 16405:0.802727 cpu time (sec):0.712119 mean cpu time (sec):0.712119
と5,6秒程度で計算終了。ウリィィィィって言いたくなるほど圧倒的に速い。
高速文字列解析の世界――データ圧縮・全文検索・テキストマイニング (確率と情報の科学)