pychembldbでつくるChEMBLウェブサービス

Flaskpychembldbを使えばChEMBLウェブサービスみたいなのは簡単に作れるよと、朝の30分くらいを使ってちょっとやってみた。

pychembldbはSQLAlchemyのラッパーなので、Flaskのほうではルーティングを設定して、ハンドラ関数用意すればいいだけ。特にFlaskはJSON化する関数が用意されているのでJSONで返すのはラク。

@app.route("/chemblws/compounds/<chembl_id>")
def compound_by_ChEMBLID(chembl_id):
    compound = chembldb.query(Molecule).filter_by(chembl_id=chembl_id).one()
    result = {...}
    return jsonify(result)

という感じでDictionaryを用意してxmlかjsonに変換して返せばいいので、とりあえずChEMBLIDを与えると対応する化合物情報を返すAPIを実装してExamplesに用意してみた。

自前でサービスを用意することのメリットは

  • 外部に情報が流れない
  • レスポンスが速い
  • 沢山投げても怒られない

ということの他に

  • 自分たちの用途に合わせて 拡張できる
  • データベースのスキーマをきちんと理解できる

という部分もあるかなと思います。例えばChEMBLウェブサービスにはジャーナルのdoiから構造リストを返すというAPIは存在しないけど、project毎にジャーナルをまとめていたりするときにはそういうAPIが用意されていると便利かもしれませんよね?

最初、ウェブサービスが返す情報は固定なのかなと思い、決め打ちで用意したのだけど、CHEMBL1CHEMBL2で返ってくるjsonのキーが違うので、valueが存在するのものをすべて返しているのかな。

もう少しちゃんと出来たらきちんとテストを書こう。

CHEMBL1

{
    "compound": {
        "acdLogd": 7.67, 
        "acdLogp": 7.67, 
        "alogp": 3.63, 
        "chemblId": "CHEMBL1", 
        "knownDrug": "No", 
        "medChemFriendly": "Yes", 
        "molecularFormula": "C32H32O8", 
        "molecularWeight": 544.59, 
        "numRo5Violations": 1, 
        "passesRuleOfThree": "No", 
        "rotatableBonds": 2, 
        "smiles": "COc1ccc2[C@@H]3[C@H](COc2c1)C(C)(C)OC4=C3C(=O)C(=O)C5=C4OC(C)(C)[C@@H]6COc7cc(OC)ccc7[C@H]56", 
        "stdInChiKey": "GHBOEFUAGSHXPO-XZOTUCIWSA-N"
    }
}

CHEMBL2

{
    "compound": {
        "acdBasicPka": 6.52, 
        "acdLogd": 2.09, 
        "acdLogp": 2.14, 
        "alogp": 2.11, 
        "chemblId": "CHEMBL2", 
        "knownDrug": "Yes", 
        "medChemFriendly": "Yes", 
        "molecularFormula": "C19H21N5O4", 
        "molecularWeight": 383.4, 
        "numRo5Violations": 0, 
        "passesRuleOfThree": "No", 
        "preferredCompoundName": "PRAZOSIN", 
        "rotatableBonds": 4, 
        "smiles": "COc1cc2nc(nc(N)c2cc1OC)N3CCN(CC3)C(=O)c4occc4", 
        "species": "NEUTRAL", 
        "stdInChiKey": "IENZQIKPVFGBNW-UHFFFAOYSA-N", 
        "synonyms": "CP-12299,Minipress,Minizide,PRAZOSIN,Prazosin"
    }
}

DMTAがいいのかそれともTADMがいいのか?

PDCAに対応するものとして創薬の文脈の中ではDMTA(Design-Make-Test-Analysis)が語られるわけだが、昨日ケミストの方にTADMでサイクルを回すのはどう思うか?と聞かれ、直感的に「それは逃げなんじゃないかなぁ、結局ランダムスクリーニングを肯定することになるじゃん」と否定的に答えてしまったのだけど、後から考えてみると根拠があまりないわなと、もう少し深く考えてみた。

DMTA(PDCA)は仮説駆動サイクルですね、これはまぁ自明。じゃぁ、CAPD(TADM)っていうのはなんなんだ?って考えてみるにこれはテスト駆動サイクルかなと思っていたのだけど、CAPDサイクルでググッてみるに、どちらかと言うと探索駆動サイクルの意味合いが強そうだ。昨日はすっかり忘れていたが、これはちょっと前に考えたことがあった

創薬系においては、そもそもテストにあたるもの(アッセイ系)を構築するのに非常にコストがかかる。同様にDoに対応するMakeにも人件費が大いにかかるため、ランダムな合成というのはリスキーだ。さらに探索という側面がある。そのため、仮説なしに合成する(これをケミストは情報取りのための合成としばしば呼ぶわけだが)、ランダムな合成をして有意義な情報が取れたことが経験上ほとんどない。

という事情もあって、昨日の質問に(半ば短絡的に)否定的な返答をしたのだと思う。

でも実際にはlead optimizationはまず間違いなくDMTAをまわすべきだけど、lead generation (finding)にはTADMサイクルをまわすべきなのだと思う。そのためには探索手法(なにを探索したいのかそのためにはどういう測定系が必要か?さらにはどういう化合物セットを揃える必要があるか?)に精通した人間(またはチーム)が必要なんだろうなぁと。残念なことに僕はそういう探索特化型のメディシナルケミストと一緒に仕事をしたことがないという背景が否定的な結論を導いたのかなと思っている。

そういうヒトと一緒に仕事を出来ればハッピーだろうなぁと思いますね。

というわけで色々考えていたら、DMTAサイクルをまわす能力よりもTADMをまわす能力のほうがずっとレアかなと。昨日の主張はちょっと短絡的だったなと反省した。

長距離索敵陣形

ProductName 進撃の巨人(5) (講談社コミックス)
諫山 創
講談社 / 450円 ( 2011-08-09 )


gWTで何をやりたいのかっていう発表をしてきた。

今週、来週と発表が続いでいるのだけど、今日は(自分の中で)アツいgWTを製薬コミュニティの中で紹介してきて、そういう手法を使って何をやりたいのかっていう話をしてきたが有意義なレスポンスをいただけて非常に良かった。特にアルゴリズム系に興味あるヒトと繋がれたのと、同じような仕事のやり方をしている人とも話が出来たので、ユーザー会に出た価値があったかなと。

というより、初めて挨拶させていただいた複数の方から「ブログ見てるで」ということを言われて話を合わせやすかったので、ブログいいよなと再認識した。

それから、WizePairZの論文のAuthorにおもろかったでと言われたのだけど、そのときは著者だと認識できてなくて後からその方が発表中にAuthorだと気づいて、「おー!論文いつも楽しく読んでますよー」と念をとばしてみた(が当方変化型なので放出レベルは低いかも)。ちなみに来週はMMP絡みの話をするというか、activity cliffをプロジェクトのインシデントしてダブラズモラサズMECE精神で機械に自動的に拾わせるという話をするのだけど、また有意義な議論が出来ればいいなぁと。

python-pptxで表をつくってみた

python-pptxがテーブル対応していたので早速使ってみた。これはヤヴァイ、スクリプトで自動化したら快適になりそうな予感がする。

が、いまのとこ6行以上を指定するとファイルが出力できない模様。

test_table

from pptx import Presentation
from pptx.util import Inches
from pychembldb import *

prs = Presentation()
title_only_slidelayout = prs.slidelayouts[5]
slide = prs.slides.add_slide(title_only_slidelayout)
shapes = slide.shapes

shapes.title.text = "10.1016/S0960-894X(01)80693-4"

rows = 6
cols = 6
left = Inches(0.5)
top = Inches(1.5)
width = Inches(8.0)
height = Inches(0.8)

tbl = shapes.add_table(rows, cols, left, top, width, height)

tbl.columns[0].width = Inches(2.0)

tbl.cell(0, 0).text = 'inchi_key'
tbl.cell(0, 1).text = 'activity'
tbl.cell(0, 2).text = 'HBD'
tbl.cell(0, 3).text = 'HBA'
tbl.cell(0, 4).text = 'PSA'
tbl.cell(0, 5).text = 'MWT'

for journal in chembldb.query(Doc).filter_by(doi="10.1016/S0960-894X(01)80693-4"):
    for assay in journal.assays[:1]:
        for i, activity in enumerate(assay.activities[:5], 1):
            tbl.cell(i, 0).text = activity.compound.molecule.structure.standard_inchi_key
            tbl.cell(i, 1).text = str(activity.published_value)
            tbl.cell(i, 2).text = str(activity.compound.molecule.property.hbd)
            tbl.cell(i, 3).text = str(activity.compound.molecule.property.hba)
            tbl.cell(i, 4).text = str(activity.compound.molecule.property.psa)
            tbl.cell(i, 5).text = str(activity.compound.molecule.property.mw_freebase)

prs.save('test.pptx')

追記 130605

というわけでSHizuoka.py #2も楽しいよと、ちょいと宣伝しておきます。

chembl_15をサポートしてバージョンもあげた(pychembldb)

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

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)。

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

ChEMBLデータベースのPythonインターフェースを書いた

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

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

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になった

化合物の入力にjava製のプラグインはもう流行んないだろうなぁと思った

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

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

あとで試す。

Graph-indexing wavelet tree (gWT)のグラフ検索が超速い

高速文字列解析の流れで、タイミングよく統数研チャンネルの資料が公開されていて、最後のほうのスライド(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はすごいなぁと思った。

OS10.8.2にopenbabel2.3.2をインストール

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

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

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