先週バージョンが上がってスキーマもかなり変更されたchembl_15に対応させました。

• pychembldb

ChEMBLはバージョンアップでおもに結合サイトまわりが強化されているみたいです。binding_siteとかpredicted_binding_domainみたいなデータが入っているのはドッキングシミュレーションを強く意識した作りになっているんでしょう。とはいえ、いい感じのOSSなドッキングシミュレーションソフトウェアがないから、遊びにくいかも。

スキーマに関して言えば、多対多はすっきりとした感じに変更されてたのでmapper書きやすくなっててよかったが、ちょいとハマったところをメモっておく

• target_dictionaryとpredicted_binding_domainsがone-to-manyになってなくてactivityに紐付いてる
• information_source doesn't exist in products table
• USAN_STEMがautoでマップできないのでとりあえずコメントしてある

あとはテストケース書くのがめんどくさかったので自動生成するようにした。今回はクラスもある程度自動生成した(utils)。

• pychembldb - Github

doctest書かないといけないなぁと思った。

ちょっと欲しくなったので探してみたらpychemblを見つけた。

しかしメンテされてないようなのでフォークして更新しようかなと思い、ソースを眺めたら面倒くさそうな感じだったので、昨日一日使って新しいのを作ってみた。

• pychembldb

declarativeでautoloadしているのでスキーマが変わってもメンテしやすいかなと思う。コーディングのほとんどの時間をテスト書くのに使ったので結構疲れた。それから、メタプログラミングはどうかなと思ったのでdbのuriは環境変数CHEMBL_URIで設定するようにしておいた(デフォルトはchemblのmysqlのインストールそのまま)。

chembl = ChEMBLDB(uri)とかやりたい場合にはFlask-SQLAlchemyのソースが参考になると思う。

ターゲットから活性値とInChiKeyをフェッチする

ターゲット名に含まれるアッセイ系から活性と対応する化合物の構造をとってくる。

from pychembldb import *
for target in chembldb.query(Target).filter_by(pref_name="Tyrosine-protein kinase ABL"):
    for assay in target.assays:
        for activity in assay.activities:
            print activity.published_value, activity.molecule.structure.standard_inchi_key

ジャーナルから活性値とInChiKeyをフェッチする

ジャーナル(doi)に含まれるアッセイ系から活性と対応する化合物の構造をとってくる。

from pychembldb import *
for journal in chembldb.query(Doc).filter_by(doi="10.1016/S0960-894X(01)80693-4"):
    for assay in journal.assays:
        for activity in assay.activities:
            print activity.published_value, activity.molecule.structure.standard_inchi_key

その他

ChEMBLのスキーマみれば大体わかるだろうし、その通りに動くでしょう。このスライドも参考になると思う。

余談ですが、Rubyだとbio-chemblがあるそうです。

今後やる予定のこと

• とりあえずテーブル名をそのままクラス名にしたせいでlookupとかDictionaryとかついててダサいのでそこら辺の名前は変える（なのでクラス名とかメソッド名はまだ変わります）
• relationが不完全なので追加していく
• one-to-oneみたいなのもone-to-manyとしてrelation張っているせいでリストになっている部分は修正する(上だとmolecule.structures[0]とかださい)
• pypiに登録する(pipでインストール出来ないと面倒くさい)

130128 追記

• one-to-oneにしたのでmolecules[0]がmolecule
• Dictionaryを除いたのでTarget,Moleculeになった

NARのDatabase IssueにMMPのサービスが載るんだとか思いつつ、SwissBioisostereにアクセスしたら、構造描画にjava製のプラグインを要求するので使うのをやめた。

• ChemWriter
• ChemDoodle

なんかが良いかなと思うが。chemdoodleのチュートリアルは結構面白いしbackbone.jsなんかと組み合わせられるんじゃないかなぁ。

あとで試す。

高速文字列解析の流れで、タイミングよく統数研チャンネルの資料が公開されていて、最後のほうのスライド(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秒程度で計算終了。ウリィィィィって言いたくなるほど圧倒的に速い。

bayonも感激したけど、それ以上にgwtはすごいなぁと思った。

高速文字列解析の世界――データ圧縮・全文検索・テキストマイニング (確率と情報の科学)
岡野原 大輔
岩波書店 / 3150円 ( 2012-12-27 )

2.3.1のままだったので、昔のログを見ながらインストール。

が、古いswig(/usr/binにあるやつ)を見ていてエラーになるので、SWIGの場所をcmakeに知らせたいのだがオプションがわからない。

調べたトコロ

cmake -L

で一覧が出るらしい。必要なのはSWIG_EXECUTABLEなので

cmake ../openbabel-2.3.1 -DPYTHON_BINDINGS=ON \
-DPYTHON_LIBRARY=/usr/local/lib/libpython2.7.dylib\
-DPYTHON_EXECUTABLE=/usr/local/bin/python \
-DEIGEN2_INCLUDE_DIR=/Users/kzfm/openbabel/eigen-eigen-2.0.17 \
-DRUN_SWIG=ON \
-DSWIG_EXECUTABLE=/usr/loca/bin/swig

でOK

見逃してた。これは超欲しい。

Bioisosteres in Medicinal Chemistry (Methods and Principles in Medicinal Chemistry)

Wiley-VCH / 15668円 ( 2012-08-22 )

opensourceでRXNかSMIRKSのあつかえるchemoinformaticsのライブラリが欲しいなぁ。

僕のフィールドには一応Bioinformatics,Chemoinformaticsも含まれているので、文字列だけでどこまでいけるのかは非常に興味がある。まぁ物理法則を無視できなくなってくると文字列処理ではどうしようもないんだけどね。

高速文字列解析の世界――データ圧縮・全文検索・テキストマイニング (確率と情報の科学)
岡野原 大輔
岩波書店 / 3150円 ( 2012-12-27 )

3章のBWTまで読んだけど、3-4のBWTの性質と復元のところがさらっと流されていてさっぱり分からなかったので、検索したらわかりやすい説明を見つけた。

• Burrows Wheeler TransformとLF mapping

が、全体的には丁寧に解説されていて、分かりやすいと思うし、読んでて楽しい。

個人的にはSMILESみたいな構造情報を文字列にしたものをBWTで扱えないかなぁと考えている(chemoinformaticsに応用してみたい)。構造変換ルールにもイディオムみたいなのあるしね。MMPを高速に検索できても嬉しいだろうし、なんか使い道がありそうな気がするんだけど。

ローカルに構築したChEMBLにアクセスする場合には生のSQLじゃなくてSQLALchemyを使ったほうが楽だ。その場合にはMetaDataのrelflectをTrueで。

from sqlalchemy import create_engine, MetaData

db = create_engine('mysql://user:password@localhost/chembl_14?charset=utf8&use_unicode=0')
metadata = MetaData(bind=db, reflect=True)

table = metadata.tables['target_dictionary']
stmt = table.select(table.c['organism'] == 'Homo sapiens')

for r in stmt.execute():
    print r.chembl_id, r.pref_name

ChEMBLって、生物種を正規化してないのね。パフォーマンスの問題かな？

Bioinformatics Programming Using Python
Mitchell L. Model
Oreilly & Associates Inc / 5075円 ( 2009-12-23 )

システム構築してるはずが、ここのところ分析多めというか120%

Canopy Clustering

ちょっと調べ物をしていてCanopy Clusteringという手法がmahoutに実装されていて面白そうだったので、論文読んでみた。

最初、自動で分割数選んでくれるのかーと勘違いしてたのだけど、荒くクラスタリングしてその後K-meansとかGreedy Agglomerative Clusteringに渡すだけらしい。分割数も荒いクラスタリングの中でAICとかBICを使って普通に求めるっぽい。

Mahoutイン・アクション
Sean Owen
オライリージャパン / 4200円 ( 2012-10-26 )

調べていく過程でIIR読めというような啓示があったが、積んであるのでそろそろ読もうかと思った。

情報検索の基礎
Christopher D.Manning
共立出版 / 8400円 ( 2012-06-23 )

Deep Learning

NNがアツいらしい。活性予測なんかでも好成績をおさめているらしいのでちょっと興味がわいた。

Merck Molecular Activity Challengeはなかなかいい試みだなぁ。日本でもどこかがやらないかなぁ(リーディングカンパニーとか)

特に興味を惹かれたのは

いくつかの分野、特に人の認識に関係するような画像、音声、言語といった分野では既存の技術を大きく超えている感があります。

ドラッグライクネスというモヤモヤ感満載の指標がRo5から置き換わるのであれば素晴らしいことかもしれない。

redmineはGitのようなVCSをセットにして使うのが当たり前なので、創薬系でもそのようなバージョン管理システムがあればいいなぁと思っている。

現状の化合物管理データベースは履歴管理をしないし、コミットログに対応するものも記録できないので、単なる構造データ管理としてしか機能しないが、化合物管理データベースのフロントにwikiかそれをモディファイしてコミットログとか履歴管理に対応する機能をもたせたものを用意すればやりたいことができそうだ。さらにAPIを用意してcliのラッパーとかpython,rubyバインディングを用意すればredmineと連携させやすいだろう。

logでログとかコミットツリーが出て

$ chemit log
commit C0000123
Author: ChemistA <XXXXX>
Date:   Sun Aug 19 10:05:53 2012 +0900

    improve solubility

commit C0000122
Author: ChemistB <YYYYY>
Date:   Sun Aug 19 10:04:58 2012 +0900

    reduce clogp  fixes #24

commit C0000121
Author: ChemistA <XXXXX>
Date:   Sat Feb 4 09:02:58 2012 +0900

代謝安定性を改善 fixes #23

それから

$ chemit diff C0000121 C000123

と打つと構造の差分とMCSが表示されようなものを作ってredmineと連携できれば素敵じゃないかなと思ったので、実証用のコードを書いてみようかな。

Redmineによるタスクマネジメント実践技法
小川 明彦
翔泳社 / 3444円 ( 2010-10-13 )