創薬アドベントカレンダー 2018

6日目の記事 #souyakuAC2018 となります

エモい＋ケモインフォ＝ケモい

— ホタペン (@hinaichigo) 2018年11月20日

それはつまり、 Python+RDKitはエモくないから、もっとエモい言語でケモインフォマティクスをやりたいということですな？

よし、やろう！

ところで、今どきのケモインフォマティクスといえばRDKit+Pythonが一般的だろうけど、その前はどうだったのかということを、自分の記憶とかブログのメモを頼りに歴史を振り返ってみようかと思う。ここ15,6年くらいの話になると思う。

LLでケモインフォ

LLをつかってケモインフォマティクスをやりたいという話が出始めたのが2000年過ぎくらいからだったと思う。バイオインフォマティクス分野ではオープンソースのツールをつかって研究をするのが一般的だったので、その流れでケモインフォマティクスにもオープンソースでやりたいというニーズが高まったように記憶している。その当時はケモインフォっていうかQSARのほうがポピュラーだったので統計解析パッケージのRとかでなんとかしようとしてたような。

ケモインフォマティックスが出たのが2005年だからそれより前は洋書しかなかった気がする（でも少なかった）

ケモインフォマティックス―予測と設計のための化学情報学
J.Gasteiger
丸善 / 18900円 ( 2005-02 )

Rでケモインフォ (2003-2007)

前の会社に入った頃はSBDDとバイオインフォ、ケモインフォを並行してやっていて、その当時はbioinformaticsといえばBioPerl一択でperlばっかり書いてた。ただQSARはRを使っていて、SVMとかMLRとかしていた気がする。大体2003-6年くらいですね。

この当時は構造描画したり、記述子を発生する唯一のオープンソースライブラリがCDKだった気がするのでみなさん他言語からCDKを呼び出そうしてた気がする。RCDKもこの当時にはあった。

結論: Rエモくない

Perlでケモインフォ (2006-2008)

Perlは学生の時から書いていてい、もちろんその当時もperlユーザーであった私はもちろんperlからCDKを呼び出すためにInline::Javaをつかうということを思いつき実践していた。

一方でPerlMolっていうPerlで構造のマニピュレーションするためのモジュールが出てきて自分の中で高まりをみせた感があって気がする。日本人で他に使っていた人がいたかは知らんけど、RECAP書いたり、人工無能作れたりとなかなかよかった気がした。

結論: perlは若干エモかったがいまさら書きたいとはあまり思わない

Ruby でケモインフォ (2006)

エモいLLといえばRuby,そしてBioRubyの流れからのChemRubyというのがあって、描画がイケていたので使おうとしたけど、メンテナであった金久研の方が忙しくなって更新されなくなってしまったので、そういうものがありましたということで残しておく

ちなみにRuby関連で読んでおく本は「メタプログラミングRuby」と「Rubyで作る奇妙なプログラミング言語」ですね。

メタプログラミングRuby 第2版
Paolo Perrotta
オライリージャパン / 3240円 ( 2015-10-10 )

Rubyで作る奇妙なプログラミング言語 ~Esoteric Language~
原 悠
毎日コミュニケーションズ / ?円 ( 2008-12-20 )

結論: Rubyのエモさがchemoinformaticsには注入されなかった

Jython + CDKでケモインフォ (2008)

私がPythonを真面目に書くようになったのが2006年ごろだったと思う。当時GAMESS FMO inputの生成スクリプトにbioperlのpdb parserをつかっていたのだけどやたら遅くて、ふとbiopythonのparserつかったらすごく速いことに感動してそのまま宗旨変えした。openbabelのbindingはperlもpythonもどっちもあったので両方つかっていた記憶がある。

Jythonだとpythonとjavaがうまく連携できるのでCDKとの相性が非常によろしかった。

けどJavaで書くのが面倒くさかったのとJavaのコード書いてもあまりテンションが上がらなかったのに加え、具体的なプロダクトにまで進みそうな仕事がなかったので、楽しんで終わっただけだった。

その他Javaで書かれたやつ(Scala, Clojure) (2009-)

ScalaやClojureを使う人が増えたり書籍も色々と出た結果、Javaで作られたchemoinformaticsのtoolkit(ChemAxonね)をScalaやClojureで取り扱う系のblogの記事とかはちょいちょいあって、「ほー、なるほどー」とか言って参考にしていたと思う。そのあたりのブログを探してみたけど今はなくなっているようだった。

• cdk-scala-examples

結論: Scala面白いですよね、マルチパラダイム言語ってのがエモい。ChemAxonの製品つかっていたら今頃Scalarになってたかも

Scalaスケーラブルプログラミング第3版
Martin Odersky
インプレス / 4968円 ( 2016-09-20 )

関数型言語でChemoinformatics (2010-?)

2006年頃から関数型言語に傾倒しだして、趣味のプログラミングとしてHaskellを触りはじめlisp,schemeに傾倒していて関数型言語でchemoinformaticsやりたくなった。

Pythonしか書かんよ？っていう人でも、とりあえずプログラミングHaskellは読んでおくとよいでしょう。モナドは出てこないので、 関数型言語おもしろい！ 感は味わえます。あとこれを読むとPythonのプログラミングがちょっと変わります。リスト内包表記とかzip系を多用する病気にかかります、あとラムダ式ね。

プログラミングHaskell
Graham Hutton
オーム社 / 2940円 ( 2009-11-11 )

次にすごいHaskell読めばいいと思います。読もう、読むべき

すごいHaskellたのしく学ぼう!
Miran Lipovača
オーム社 / 2940円 ( 2012-05-23 )

で、このあたりが分かりだすと色々やりたくなるわけですね。以下の論文読むと色々書いてあるけど、型安全とか参照透過性とか遅延評価バンザイとか実際コード書いていると欲しくなる。

• chemf: A purely functional chemistry toolkit

で、Ouchなどを追いかけていたけどなかなか仕事で使う感じにはならなかったので結局python + openbabel(仕事ではOE toolkit)でコードを書いてた。

結論: Haskellはエモいけど、使えるライブラリがなかった ;-)

Pythonでケモインフォマティクス (2008-)

そんなわけで2008頃からPythonでケモインフォマティクスのコードを書くようになったんだけど、その当時はPythonistaも今ほど多くはなかったですね。暇だったのでこのブログとは別にケモインフォクックブック書いてた。ちょうどこの時期にopenbabelのpythonラッパーであるpybelが構造描画に対応したのでオープンソースでも使えるようになってきたと記憶している。でもscikit-learnはなかったのでpythonからRにアクセスできるRpyつかっていた。その当時は論文の実装Rだったし、PythonからR呼び出すのが都合が良かったと思う。scikit-learnは割と最近なんじゃないですかね。

それからRDKitがどのくらいから盛り上がってきたのかを調べるには日本のパイオニアとして作者にも認識されているpen先生のブログを追えばいいのですが、それによると大体2012年くらいですね。ただこの頃はopenbabel派が多かったのでRDKitが主流になるのはもう少し後のはずです。

一応openbabelとRDKitの違いを述べておくと、openbabel(OpenEyeも含む)はケミストリー寄りで、RDKitはケモインフォマティクス寄りです。なにが違うのと言われると、「違うよ、ぜんぜん違うよ！！！！」くらい違うんですが、openbabelは量子化学計算の結果を格納することを前提とした作りになっているのに対し、RDKitはそこは考慮されてないんですよね。

結論: Pythonエモくない、固い、カッチカチ

エモいケモインフォマティクス

ここから本題、前置き長かったわ。

lispでケモインフォやれば完璧 https://t.co/NFBKDKEImS RT @hinaichigo: エモい＋ケモインフォ＝ケモい

— kzfm (@fmkz___) 2018年11月20日

エモいプログラミング言語といえば関数型でhaskellとかlispですが実はpythonによるlisp実装の一つにhyというものがあります。

インストールも簡単、みんな大好きpipで

pip install git+https://github.com/hylang/hy.git

と叩くだけでOK

入門用の参考リンクを3つほどあげておきます。

• Pythonで書かれたイカしたLisp: Hy
• ゼロから始めるHy
• python使いのためのhy基本文法まとめ

Pythonで書かれた活性予測のコードも

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import DataStructs
import numpy as np
from sklearn.metrics import r2_score
from math import sqrt
from sklearn.svm import SVR
from sklearn.model_selection import train_test_split

def svm(x_train, x_test, y_train, y_test):
    clf = SVR(kernel='linear').fit(x_train, y_train)
    y_pred = clf.predict(x_test)
    r2 = r2_score(y_test, y_pred)
    print('SVM: R2: {0:f}, RMSE:{1:f}'.format(r2, rmse))

def sdf_to_desc(sdf_file):
    fps = []
    targets = []
    nfps = []

    for mol in Chem.SDMolSupplier(sdf_file):
        fps.append(AllChem.GetMorganFingerprintAsBitVect(mol, 2))  # fingerprint
        targets.append(float(mol.GetProp("ACTIVITY")))  # activity

    for fp in fps:
        nfp = np.zeros((1,))
        DataStructs.ConvertToNumpyArray(fp, nfp)
        nfps.append(nfp)

    return (np.array(nfps), np.array(targets))

if __name__ == '__main__':
    descs, targets = sdf_to_desc("CHEMBL1827733_5HT2A.sdf")
    x_train, x_test, y_train, y_test = cross_validation.train_test_split(
        fps,
        acts,
        test_size=0.1)
    svm(x_train, x_test, y_train, y_test)

hyで書きなおすとあら不思議、括弧の多さがエモさの証、こんなにエモくなりました!

(import [rdkit [Chem]])
(import [rdkit.Chem [AllChem]])
(import [rdkit.Chem [DataStructs]])
(import [math [log10]])
(import [math [sqrt]])
(import [numpy :as np])
(import [sklearn.model_selection [train_test_split]])
(import [sklearn.svm [SVR]])
(import [sklearn.metrics [r2_score]])

(setv sdf_file "CHEMBL1827733_5HT2A.sdf")
(setv acts (map (fn [mol] (- 9 (log10 (float (mol.GetProp "ACTIVITY")))))
    (Chem.SDMolSupplier sdf_file)))
(setv fps (map (fn [mol] (AllChem.GetMorganFingerprintAsBitVect mol 2))
    (Chem.SDMolSupplier sdf_file)))
(setv nfps [])
(for [fp fps] (setv nfp (np.zeros '(1, )))
    (DataStructs.ConvertToNumpyArray fp nfp)
    (nfps.append nfp))
(setv Y (np.array (list acts)))
(setv X (np.array nfps))
(setv D (train_test_split X Y :test_size 0.1)) ;D:x_train, x_test, y_train, y_test
(setv svr (SVR))
(setv clf (svr.fit (get D 0) (get D 2)))
(setv y_pred (clf.predict (get D 1)))
(print "r2:" (r2_score (get D 3) y_pred))

マクロも使えるみたいだし、積んであったOn Lisp読むモチベーションが高まったわ。

On Lisp
ポール グレアム
オーム社 / 3990円 ( 2007-03 )

結論：AI創薬をうたうならば括弧を恐れず本物の人工知能言語であるlispをつかってエモいケモインフォをやるべき。