創薬アドベントカレンダー 2018
6日目の記事 #souyakuAC2018 となります
それはつまり、 Python+RDKitはエモくないから、もっとエモい言語でケモインフォマティクスをやりたいということですな？
よし、やろう！
ところで、今どきのケモインフォマティクスといえばRDKit+Pythonが一般的だろうけど、その前はどうだったのかということを、自分の記憶とかブログのメモを頼りに歴史を振り返ってみようかと思う。ここ15,6年くらいの話になると思う。
LLでケモインフォ
LLをつかってケモインフォマティクスをやりたいという話が出始めたのが2000年過ぎくらいからだったと思う。バイオインフォマティクス分野ではオープンソースのツールをつかって研究をするのが一般的だったので、その流れでケモインフォマティクスにもオープンソースでやりたいというニーズが高まったように記憶している。その当時はケモインフォっていうかQSARのほうがポピュラーだったので統計解析パッケージのRとかでなんとかしようとしてたような。
ケモインフォマティックスが出たのが2005年だからそれより前は洋書しかなかった気がする（でも少なかった）
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のエモさが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の記事とかはちょいちょいあって、「ほー、なるほどー」とか言って参考にしていたと思う。そのあたりのブログを探してみたけど今はなくなっているようだった。
結論: Scala面白いですよね、マルチパラダイム言語ってのがエモい。ChemAxonの製品つかっていたら今頃Scalarになってたかも
関数型言語でChemoinformatics (2010-?)
2006年頃から関数型言語に傾倒しだして、趣味のプログラミングとしてHaskellを触りはじめlisp,schemeに傾倒していて関数型言語でchemoinformaticsやりたくなった。
Pythonしか書かんよ？っていう人でも、とりあえずプログラミングHaskellは読んでおくとよいでしょう。モナドは出てこないので、 関数型言語おもしろい！ 感は味わえます。あとこれを読むとPythonのプログラミングがちょっと変わります。リスト内包表記とかzip系を多用する病気にかかります、あとラムダ式ね。
次にすごいHaskell読めばいいと思います。読もう、読むべき
で、このあたりが分かりだすと色々やりたくなるわけですね。以下の論文読むと色々書いてあるけど、型安全とか参照透過性とか遅延評価バンザイとか実際コード書いていると欲しくなる。
で、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エモくない、固い、カッチカチ
エモいケモインフォマティクス
ここから本題、前置き長かったわ。
エモいプログラミング言語といえば関数型でhaskellとかlispですが実はpythonによるlisp実装の一つにhyというものがあります。
インストールも簡単、みんな大好きpipで
pip install git+https://github.com/hylang/hy.git
と叩くだけでOK
入門用の参考リンクを3つほどあげておきます。
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をつかってエモいケモインフォをやるべき。
ビジョナリー・カンパニー２ 飛躍の法則
Waystone
The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems
ピューと吹く！ジャガー モノクロ版
ケモインフォマティックス―予測と設計のための化学情報学
メタプログラミングRuby 第2版
Rubyで作る奇妙なプログラミング言語 ~Esoteric Language~
Scalaスケーラブルプログラミング第3版
プログラミングHaskell
すごいHaskellたのしく学ぼう!