そういえば多数の化合物からなるsdfからExactMassってどうやって計算したらいいのかね？と尋ねられて、

うーんそれはね、perlでChemistry::OpenBabelモジュールを使ってだなぁ、ごにょごにょとやればいいんだよ、ほら楽勝だろ。

とか言ってみたんだけど、別にpythonでもいいので書いてみた。

import openbabel as ob

obconv = ob.OBConversion()
obconv.SetInFormat("sdf");
obmol = ob.OBMol()

notatend = obconv.ReadFile(obmol, "test.sdf");
while notatend:
    print obmol.GetExactMass()
    obmol = ob.OBMol()
    notatend = obconv.Read(obmol)

たぶんPlutoだったら

from pluto import *

for mol in Molecule.readfile("test.sdf"):
    print mol.ExactMass

と書けるようになるはず。

でも今夜はOBChemTsfmとかいうクラスがねーとかいう謎の深みにはまりかけ、SWIGやらC++にやられっぱなしでそれどころじゃないのであった。

こんな感じで。 titleとかはmol.title("newtitle")にしたほうがいいのかこのままでいいのか決めかねているのでそのうちちゃんと考える。

エタンをエチレンに変えてみた

>>> from pluto import *
>>> mol = Molecule.readstring("CC\tethane")
>>> mol.bond([1,2]).order(2)
>>> mol.title = "ethylene"
>>> mol.writestring()
'C=C\tethylene\n'

あとはmolecule同士のの結合とatomの追加と削除を実装すればよさそう

PythonでもPerlMolみたいに化学反応を手軽に扱えるツールが欲しいなぁという欲求があったので書いてみることにした。PybelみたいなOpenBabelのラッパーだけど。

そういえば、PyMOLやPyChemという、それっぽい名前は既に使われてるんですよな。

で、アトムにちなんでみた。

PLUTO (1) (ビッグコミックス)
浦沢 直樹
小学館 / ￥ 550 ()
在庫あり。

でも、国内しかわからなんだろうな。海外だとastro boyなんでしょ？まぁいいや。

>>> from pluto import *
>>> mol = Molecule.readstring("c1ccccc1O")
>>> mol.title = "phenol"
>>> mol.writestring()
'c1c(cccc1)O\tphenol\n'

とりあえず、readとwriteまではできた。SMARTSのクラスも作ってあるのだけどアトムばっか返ってきてボンドの返し方がよく分からん。

diels-alder反応ぐらいまで書けるようにするとこまではとっととやる。

なかなか面白い内容だった。

Pybel: a Python wrapper for the OpenBabel cheminformatics toolkit

Pybelボンドとか扱えないような感じなんで、ちょっと複雑なことやろうとすると結局Openbabelのインターフェースからアクセスしちゃうんだよなー。OpenBabelのラッパーとしてPerlMolみたいなの書けばいいんだろうけど。

とか思いながら、OpenBabelのAPIなどみてたら、PerlMolと同等なものを書けそうな気がしてきたので少しやってみることにする。

というわけで、分子量とかformulaとか出せるようにしておいた。

CCC       propane
CC        ethane
c1ccccc1  benzene
CO        methanol

というsmilesを

babel -ismi test.smi -osdf testsdf

で変換したファイルを使った。ExactMassを計算する場合は

>>> from pluto import *
>>> mols = Molecule.readfile("test.sdf")
>>> for mol in mols:
...   print mol.exactmass
...
44.062600256
30.046950192
78.046950192
32.026214748

ついでに分子量と化学式も

>>> from pluto import *
>>> mols = Molecule.readfile("test.sdf")
>>> for mol in mols:
...   print mol.formula, mol.molwt
...
C3H8 44.09562
C2H6 30.06904
C6H6 78.11184
CH4O 32.04186

openbabelにはPerlMolでいうところのcombineみたいなメソッドがないので反応させるところで悩み中。

BKChemという分子構造式エディタがあって、pythonで書かれているうえに、スクリーンショットみると構造なんかも綺麗。で、バッチ処理もできるらしいので、ちょいとこれを使ってみるかなとWindowsのバイナリ配布のやつを入れてみた。

だが、ありがちなことに、このバイナリ版はバッチ処理が出来なかったのでソースから入れたけどpycairoを要求してきたのでGTK関連を一式インストール。

GTK+のバイナリはgtk.orgからAll-in-one bundleってのを落としてきて、c:\gtk2とでも名前を変えてパスを切っておく。

で、python bindingはここのを使った。というよりeasy_installではうまく入らなかった。

これでBKChemはすんなりインストールできる。実行ファイルはC:\Python25\Lib\site-packages\bkchem\bkchem.pyだった。

bkchem.py -b batch.py args

という書式でバッチ処理が行われるようになったのでsdf2pngみたいな処理を行わせればcdkの代わりに使えるはず。

pycairo自体も結構面白そうなのでsampleをやってみた。

他にもチュートリアルがあるのであとで読もう。

Fetching PDB files remotely in pure Python codeを見つけて、 数行でpdbのファイルをフェッチできるなんてやるなRCSBとpythonのコンボとか思ったんだけど、コメントに「生でフェッチすんのは環境に悪いでよ」とか書いてあった。

確かにPDB界ではgzip圧縮したpdbファイルをとってきてローカルで展開すんのが昔からのナラワシだよなと、試しにgzとかZをくっつけてwgetしてみるとその通りのファイルがダウンロードできた。をを。

で、おーこれは、数行スクリプトにgzipモジュールかませばいいだけなんじゃなかろうかと思て書き書きしてみたが、これってファイルしかとれないのかな?

zlibならどうだといじってみたけど、

>>> import urllib
>>> import zlib
>>> url = 'http://www.rcsb.org/pdb/files/1ab6.pdb.gz'
>>> text = urllib.urlopen(url).read()
>>> zlib.decompress(text)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
zlib.error: Error -3 while decompressing data: incorrect header check

なんかヘッダがーっておこられた。やっぱgzipでなんとかすべきなのか? ともうちっと追いかけたらGzipFileのほうを使えばよいらしいことに辿りついた。

>>> import gzip, StringIO, urllib
>>> url = 'http://www.rcsb.org/pdb/files/1ab6.pdb.gz'
>>> content = urllib.urlopen(url).read()
>>> sf = StringIO.StringIO(content)
>>> dec = gzip.GzipFile(fileobj=sf)
>>> data = dec.read()

ちと長いがこれで読めた。というわけで、

>>> import urllib, StringIO, gzip
>>> def fetch_pdb(id):
...     url = 'http://www.rcsb.org/pdb/files/%s.pdb.gz' % id
...     content = urllib.urlopen(url).read()
...     sf = StringIO.StringIO(content)
...     return gzip.GzipFile(fileobj=sf).read()

ちょっと長くなったけど帯域にやさしい。

なつたん: Python でSICP4.3　Nondeterministic Computingを読んだら、非決定計算が割りとすんなり理解できた。去年はambわかんねと思ったが、分かってみるとあー非決定計算ねと思えた。

あとジェネレータが便利じゃないかと。それにしても遅延評価みたいなもんかなと思ったら、似たようなシチュエーションでも使えるのね。

再帰とジェネレータ

が面白く読めた。

そういえば非決定計算を使った構造活性相関解析例なんかはあるのだろうか？prologなんかを使った例は見たことあるけど、実用例じゃなかったからな。pythonでかけるなら書いてみたいなぁ。

pylonsで対話的にモデルを操作したい場合にはdevelopment.iniのあるディレクトリで

paster shell

とすればよいらしい。modelっていう変数に色々入っているので

session = model.Session
mol = model.Mol
for m in session.query(mol).filter(mol.MW > 199):
    print m

などとやるべしナ感じ。

次の二つの章が特に面白かった。

• 5 Pythonから見たオブジェクト指向
• 7 Jythonの応用例

Jythonプログラミング
西尾 泰和
毎日コミュニケーションズ / ￥ 3,150 ()
在庫あり。

単なるjythonの使い方の本というよりは、jythonのような複数の言語を混ぜ合わせた処理系を通して見えるものを感じ取りましょう的な内容だったように思う。あとその場合に固さやわらかさの長所をうまくのばしてやるようなプログラミングのお作法はこんな感じですよみたいな指針。

perlでもInline系のモジュールを使えば直接他の言語を埋め込んだりできて他言語のライブラリを使えるけど、言語自体が両方の世界に干渉できるともっと楽ができる。

• Inline::JavaでCDKのニ次元構造立ち上げをperlから
• Statistics::R なんかいまいち

Inline::みたいなやり方でも、jythonみたいなシームレスな言語でも結局複数の言語を理解してないとうまいことやれないので、その分習得コストはかかると思うんだけど、シームレスな言語使うと楽しいのはスケッチするような感覚が得られることなのかなと思ってる。

そういう意味ではGainerも似たような感じかも。

Make: Technology on Your Time Volume 04
オライリー・ジャパン
オライリージャパン / ￥ 1,575 ()
在庫あり。

あとシームレスな言語を使うっていうことはそれにあわせてライブラリのほうも対応していかないといくのが望ましいのかなとPybel as a generic API for cheminformatics libraries - proof of concept using CDKというエントリを見てふと思った。

chemoinformaticsだとopenbabelっていうc++の大きいライブラリとCDKっていうjavaのライブラリがあって、その上にpython(とruby)がまたがっている感じなのだけど、さらにこれらでいじくったデータを解析にもっていくためにRが必要だったりとかするのでRpyが役に立ったりする。