Tag Archives: RDKit

RDKit & Scikitlearn

Let’s try to build model using RDKit and Scikit-learn . Scikit-learn is simple and efficient tools for data mining and data analysis. I referred this page. here At first, build model using molecular descriptors . Here we go. OK! Next, build model using molecular fingerprint . Let’s try ! RF method using molecular descriptors showedContinue reading “RDKit & Scikitlearn”

CHEMBL_MMP

またしてもMMPねた RDKit のMMPモジュールを使ってCHEMBLのデータを全部くっつけるバージョン。 ただし今回はkinase_sarfariのみ。 準備として https://www.ebi.ac.uk/chembl/sarfari/kinasesarfari ここからデータを持ってきます。 ks_compound.txt ks_bioactivity.txt でディレクトリにrfrag.py, index.pyと上のデータを入れます。 くっつけるのは、、 chembl_mmp.pyとして こんな感じ ターゲットごとのMMPができます。

RECAPとか

RDKitのRECAPルールを使うと 分子をある一定の結合で切断できます。 ライブラリーにアプライして得られる頻度を見ると これよく使ってるよね。みたいなのが分かるかもと思います。 こちらのブログを拝見させていただき、pytagcloudなるものを知りました。 使うのにSDFとかpygamesとか入れないとなんなくて面倒でした。 が一応環境ができたのでRECAPで分子を切断した後 HTMLに出力するようなスクリプトを書いてみました。 ちなみにほとんどは、githubにあがっているpytagcloudのunittestのコードを ベースにしています。 で実行します。 htmlのテンプレートは https://github.com/atizo/PyTagCloud/tree/master/pytagcloud/test/web のtemplate.htmlです。 で、 カレントディレクトリにwebというフォルダを作ってそこに入れておきます。 outというフォルダにtagcloudなsmilesが作成されます。多分。 分子数が多いとRECAPで結構時間を食うし、fragmentは1っこでもカウントしてるので 適当な数以上になったらタグにするとかしないと重たい感じ。 リンクになってるからマウスオーバーしたら親化合物が見えるとか、 smi2jpgみたいになると素敵ですが、、、ちょっと技術力不足、、、 JSのらいぶらりだけどTagCloudjsというの今日教えてもらったが便利そう。

雑記

これ早いらしいです。 いろんなページで見るのですがどうなんだろうって、興味はあったけど 手を出していませんでした。 が、どうもインストールはターボール落としておしまいだそうで。 ちょっとやってみたくなったので早速、、 速度だいぶ違いました。 少しいじってみたくなりました。

Lillyのフィルタ

J.Med.Chem. 2012, 55, 9763-9772 Lillyのこれまで積み重ねてきたフィルタに関しての報告です。 ソースコードがGitHubにあがっているので、ちょっと使ってみました。 でインストールはよろしくやってくれます。簡単です。 続いてデータセットの準備をします。CHEMBLからDPP4阻害剤のデータを取ってきました。 入力に使うフォーマットはsmiles tab idにしないとだめなので 今度はPANDASで加工します。 で確認します。 ということでファイルができました。 実際に計算してみましょう 結果はsmiles / id / demerit / detailsのような出力です。 17クラスからなる275のルールから構成されるフィルターで バイオレーションによってDemeritのスレッシュフォルドが決まります。 デフォルトは100でリラックス（ゆるめ）で160をカットオフにしています。 テスト用のデータが４万件くらいでしたが、手元のMBAで数十秒でした。 プロジェクトが進んでからのアプライは？？？ですが、 市販のセットやHTS用のセットに一度当てるのは面白そうですね。

RDKitで2d_pharmacophore

RDKitにある２Dファーマコフォアを使って シミラリティマトリックスを作るスクリプトを書いてみた。 　下のデフォルト設定では3000ビット以上のFPになるので 適当なbinに指定し直した方が計算が速くなるかも。 分子によってはGitBitIdxでエラーになるのでそれを回避するため try/except処理にした。エラーが出る分子については後で原因を追う。 TCの計算はお決まりの書き方でやると、エラーにうまく対処できない。 　見にくいけどdict型でindxとFPの関係を保つようにしてみた。 よくよく考えると、SVMとかするならFPをbitstringで返しておしまいにした方がいいと気がついた。 使い方 $ python p4core.py で入力SDFに関してシミラリティマトリックスを出力される。 ECFP/FCFP/MACCSのTC 諸々と合わせて出力させると違いが分かるかな。

RDKit-MMPA2

今日は昼休み中にRDKitのMMPAスクリプトをいじってました。 kzfmさんのブログを拝見させていただいたのがきっかけです。 RDKit_2013_03_2/Contrib/mmpa配下のrfag.pyの切断ルールをRECAP風にアレンジしてみました。 209行目以降 を とします。 これでエーテル、アミド、環状アミン、エステル、芳香環−芳香環、芳香族n−芳香環、スルフォンアミド、芳香環−N−芳香環（buchwald_amination）のボンドのサーチをかけてマッチした部分の結合に関わるアトムインデックスをゲッツします。 エーテルは飽和炭素ー酸素ー飽和炭素のインデックスをとるのでインデックスは前のペアと 後ろのペアを返すようにしました。 同様に芳香環−N−芳香環も同じです。 アミドやエステルはC(=O)N, C(=O)Oのインデックスになるので０と２のインデックスだけとるようにリスト内包表記にします。 で、後はカットするのですがデフォルトのトリプルカットですとバレンスがおかしくなる可能性があるのでシングルカットにしました。 二回目以降をコメントアウトしただけですけど。 シングルカットなのである意味Rグループ分解と同じと言えば同じですが 合成する側から見ると、一気に処理できるのはいいかも。 フラッグメチルやハロゲンの効果はこのスクリプトからは見えませんけど。 あ、これの後は一緒にあるindexing.pyをいじらず実行するとRECAP風MMPができます。 SMART表記が未だにしっくりこないのでその辺が課題です。 ちょっとコード自身もきれいじゃないですね。

RDKit / PostgresSQL

RDKitはPostgresSQLと連携させて構造情報を色々と扱うことができます。 ChemBL DBはPostgresSQL用のファイルが利用できますし、 ローカルにDBを作って作業すると、自分の幅が少し広がりそうです。 SQLなんかの勉強にもなる？ 取りあえずやってみました。 まずpostgresql9.2.4をビルドしました。 マニュアルに従って、、、、 でOK 続いてRDKitの設定 でエラーが出なかったら成功。 途中で何個かwarningでましたが、一応自分の環境では８個のテスト全部パスしました。 ちなみに RDKit2013_03_2 PostgreSQL 9.2.4 BOOST 1.53.0でビルドしました。 で実際にデーターベースを作ってみます。 先ほどのフォルダの直下にdataという１０００件のsmilesが入ったファイルがあったのでこれを使ってみます。 これでテスト用のDBを作りました。 続いてファイルを読み込みます。 copy xxx from yyy;でファイルを取り込みます。デフォルトの区切りはタブなのでそのまま。 続いてrdkitの機能を入れる（という表現がいいのか不明ですが） でmolsテーブルができました。 次にプロパティの計算をしてみます。 1000件の分子量とlogpの計算が２秒弱なので早いですね。（と思うんですがどうでしょう。）

KNIME-RDKIT-R

RDKitとPythonを使ってデータを加工して Rで解析するというケースを頻繁にやっていると全部つながるといいなーと思います。 RPY2を使ってうまく繋げるのもいいですがKNIMEを使うとそれぞれの連結ができて便利かも。 ということでトライしてみました。 KNIME上のRノードでは全部Rという名前のオブジェクトに対して処理をするようで それを理解するのに時間がかかりました。 例えばPLS回帰をする場合はこんな感じです。 Rノードおのおのでplsパッケージ読み込まないと、 おいおいmvrなんて扱わないぜ。とおこられます。 で作ったのはこんな感じです。 一個はモデルを作ってモデルと実測を比較するため もう一個はモデルに当てるためのフローです。 RDKITノードでディスクリプタは簡単にとれるし、そのあとのモデリングは上のようなスニペット書いとけば後はよろしくな感じでいいので、楽できそうです。 ローカルで全部やってるからパフォーマンスが悪い、、、 まあいいか。

matplotlibでPCA

web.pyを使ってサービスを提供できるようになったので 何か使えないかな〜と考えて見た。 任意のSDFがあったらPCAやって結果を返すとかあると PCAってなに？みたいな人もspotfireで図が見えるしいいかもと思った。 PCRの方がより良いと教えていただいたのでこちらの実装も考えよう。 さて、自分がやるならRを使うのだが、web.pyに入れるとなると全部pythonがいい。 rpy2入れたしこれで。と思っていたら動かないで肩すかし、、どうも社内の環境のwinで構築したR2.10はだめなようだ RDKitのモジュールはもうメンテしてないし、matplotlibがいいじゃないと言われたのでそれで取りあえず書いてみました。 化合物数＞ビット（2048）の場合はそのままでいいみたいですが、 化合物数＜ビットの場合はそのままではエラーになるので転置します。 で最後にもう一回スコアを転置して第三成分までを出力。 手ものとのデータでやったらあまりぱっとしない感じ。 もちっとデータ集めてやってみる。

RDKitでMatched Pair Analysis

Pythonで動いて多彩な機能を持っておりながらフリーである。 ということで自分はRDKitが大好きです。 次のバージョンのβ版がリリースされているのですがどうもMMPAアルゴリズムが実装されるようです。 たくさんの情報に埋もれた何かを見るにはいいかもしれないというのと、コードに興味があるのとで早速使ってみます。 SVMはここで公開されているので 待てない人は使ってくれよ！という話です。 ちなみにリファレンスはこちら。 アルゴリズム自体は前から報告されていたのですがソースコードは今回が初です。 まだソースをちゃんと読んでないですが使ってみます。 取りあえずChemblにいってABL/Humanで1504化合物のSDFをとってきました。 入力は smiles\tidのテキストなので そのようにかえます。 でファイルができたんでフラグメント化します frag.pyでsmilesの分子をフラグメントに切断します。 1504分子から37177のフラグメントができました。 元分子　ID フラグメントみたいなテキストがchemblfrag.txtにできます。 この辺をpythonで考えると辞書型で作る感じですよね。 で最後に ペアを作ります。 結局886個のペアができました。 構造1,構造2,ID1,ID2,変換のsmirkes みたいになっています。 Spotfireに突っ込めば全部構造になるので良さそうですね。 切断部分で1504化合物で5分かかっているので大規模にやるのは時間がかかりそうだ。 合成した化合物を逐次フラグメント化してDBに入れておけばマッチングは早そうだからきびきび行くかな。 web.pyでがんばってwebインタフェースができればみんながアップしたら的なものできそうだけど。 POSTで受けたファイルのハンドリングがよくわからんからweb.pyも理解せねば こういったこと、社内で言ってもきっと理解されないなぁ 何か考えるかな〜

特許をパースする-2

WIPOのpatentscopeで作成した RSSからタイトルとサマリを抜き取って 英数字以外で区切った単語リストを作るスクリプト。 sys.argv[1]にRSSのリストをしていすると。 パテント番号と単語の数を出力する。 が ターゲットや疾患がばっちりサマリに載ってる訳でもないしなあ。 ドキュメントでクラスタライングするのは難しいかな。

WIPOのRSSからの情報取得

PATENTSCOPEはクエリをくんだ結果をRSSで出力できるので 色々解析に使えるだろうと思われるのですが いかんせん修行不足の身、最近読んでる本を参考に タイトル、番号、リンク、サマリを出すスクリプトを書いてみました。 エンコードでエラー出まくりだったので全部uft_8にエンコードしたのですがこれほんとにいいのか不明。 取りあえずPCT&PDE4 @ FPをクエリにして pde4.txtに http://patentscope.wipo.int/search/rss.jsf?query=FP%3APDE4+&office=+%28OF%3Awo%29&rss=true&sortOption=Pub+Date+Desc を保存しときます。 で patentparser.pyを書いて 取りあえず動く。 リダイレクションをファイルにして ということで出力先をもう少し整形しようと思う。 プロキシハンドラも入れないとだめですね。

USRCATを実装してみる

USRCAT = Ultrafast Shape recognition with pharmacophoric constrains RDKit のMLで興味深い報告があったので 実装してみました。 OETookit以外にRDKitもオッケーというのはとても魅力的です。 ShapeベースのバーチャルスクリーニングにてScaffold Hopingにも使える可能性があることや ElectroShapeと比較しても遜色ないEFを出している点はちょっと引かれます。 まずここからダウンロードします。 あとはマニュアルに従って $ cd /usr/local/src $ mv ~/Downloads/usrcat* $ cd /usr/loca/usrcat $ sudo python setup.py install でOK PYTHONPATHに上記のフォルダを指定しておけば、何所からでも呼び出せます。 その後、 ①まずSDF読み込み ②三次元配座の発生UFFによるエネルギーの計算 ③usrcatモジュールでのモーメントの算出 ④類似性の計算 を行わせてみました。 usrcatのgenerate_momentsはmolオブジェクトが持つ全配座のモーメントを計算します。 次のsimilarity(m1,m2)で発生させたモーメント同士の類似性をみていています。 モーメントは配座の数だけあるようですが m1は最初のconfId=0のものだけしか使わないようなので 最低エネルギーのは配座を持ってくるか、PDBからリガンドを持ってくるのがいいのかと思います。 下の例はそのまま使っているので あまりいい例じゃないです。 取りあえず動いたレベル。 コードも助長だったり見にくかったり。。。 ChEMBLのデータで何か解析してみないと。

RDKitで配座を発生させてエネルギーを計算してみる

RDKit2012_09はしばらくβ版で確認していましたが正式版がリリースされました。 今回はMCSを実装した点が今までと変わった点のように思います。 残念ながらまだwin64bit用のバイナリは無いようです。取りあえず自分のMBAではビルドもできて動作も確認できました。 MCS関連はまた今度試すとして今回は、SDFから配座を発生させて、UFFによる最適化後のエネルギーを計算させるスクリプトを書きました。 AllChem.EmbedMultipleConfsで分子に対して１０個の配座を発生させます。 AllChem.UFFOptimizeMoleculeにて構造を最適化します（1000回に設定）０なら最適化が終わっていて １、がかえってくる場合は1000以上の数値を設定した方がいいと思われますが今回次のステップにはここのクライテリアを抜いてしまいました。 ff = AllChem.UFFGetMoleculeForceField(mol, confId = cid) ff.Minimize() energy=ff.CalcEnergy() にて最適化後のエネルギーを受け取ります。 今回は出力を全部にしたのですが、最低のものだけを取り出したいのです、、、 リストのminでは数値が出てくるけどインデックスがとれないのでどうしたらいいんだろう。 list.index(min_energy)で行けそうだから後で追加しよう。

align-itを使ってみた

以前からファーマコフォアベースで分子をアライメントするツールをささっとかけると 良いな〜って思ってました。 社内リソースを使おうと思うと、ライセンスの限りがあるのでオンデマンドに解析できなかったりするためです。 RDKitでやっている事例があったのでそれを参考に書いてみたのですが、途中でだんだん複雑になってしまい、ちょっとそのまま放置しています。 そこで、もう少し簡単にということでalign-itを使ってみました。 これは前にも紹介したsilicos-itから公開されているオープンソースです。 Openbabelがカバーするフォーマットを入出力に使えるので、特定のフォーマットに変換するという作業から解放されてよいと思います。 インストールはマニュアルに従って でよいと思います。ubuntuにパッケージマネージャーでopenbabel入れてある状態でやると openbabel2 無い。とエラーがでました。どうもパッケージマネージャーではな入らないみたいですのでこの場合は、自分でopenbabelを最初にビルドしてやるといいかと思います。 続いてテスト用のファイルの準備です。 今回は面倒なんでpubchemでEGFR阻害剤のデータ21化合物をとってきました。 二個目の分子をリファレンス(ref.sdf)としてアライメントします。 結果 .pharにはフィーチャー　ファーマコフォアの三次元座標　ガウシアン体積のα値　方向性があるかないかのブール値　そのベクトル のタブ区切りがかえってきます。このデータって何で見れるんだろう？？？ またオプションの -sはスコアを返すので となりそれぞれの列はマニュアルにあるように column Content —— ——————————————————————— 1 Id of the reference structure 2 Maximum volume of the reference structure 3 Id of the database structure 4 Maximum volume of the database structure 5 Maximum volume overlap ofContinue reading “align-itを使ってみた”

特許をパースする

不勉強にてxmlって ブラウザで見られるけど何これ？ と以前思ってました。、、、恥ずかしい よくよく調べるととても便利そうです。 で、xml形式の特許があると解析に使えそうと思うわけです。 集めてきた特許から、番号や、優先日、出願人とかの情報をささっと一覧にできると いいよなーと思っていたこともあったので、 google先生に聞くとそんなことずいぶん前からなされているようでした。。。 さて、xml形式の特許が結構簡単に入手できる環境が整ったので テスト的にスクリプトを作ってみました。 これは同じフォルダにある.xml形式のファイルをを全部読み込んで出願人とかをタブ区切りで出力します。 　フランス語とか、エンコでこけるのでアブスとは英語に限定しました。 相当行けてないスクリプト感ただよいますが、一応動いた。のでモチっと改良する予定 xmlの階層構造がそもそも勉強に使っていた例に比べるとかなり入り組んでいたのでまだ不具合たくさん出そうです。 がんばると、テキストなのでバラバラにして複数の出願人がどれくらいの頻度出るかとか、簡単に見れそうです。

shape-itを使ってみた

2Dの分子の類似性を比較する場合、各種FP計算を実施してあれこれ考えるわけですが、 社内的都合にて3Dの場合、ささっと使えるツールが無い訳です。世の中的にはOpenEye社のROCSが使われることが多そうな気がします。文献出の引用例や、ユーザー会のプレゼンをみてもかなり良さげです。 お気に入りのRDKitにはShapeTanimotoメソッッドがありますが、アライメントは別途かける必要があるのでちょいと面倒でございました。 そんなおり、某ワークショップで、silicos-itのshape-itは？と教えていただきました（ありがとうございます）。 winではちゃちゃっとBuildするのは無理そうなんでまずはMacでやってみます。 shape-itのマニュアルを読んでみるとベースにOpenBabelを使っているのでOpenBabelがフォローするファイルフォーマットならなんでもOKという点も良い感じです。 Buildはマニュアルに従って, で準備は完了です。 続いてテスト用のファイルの準備です。全部openbabelでやればいいという噂もありますが、こちらは自分のスキルの都合上RDKitで取りあえずのファイルを作ります。テストファイルはpubchemからもらったDYRK1Kのデータ。 2Dの分子から3D配座を発生させて取りあえずUFFで一応最適化しておきます。 全部の準備が終わったので 実行してみました。 ここでは先頭の分子をreferenceとして使います。 結果 volumeだけの比較になりますが、ちゃんとアライメントもかけて出力してくれそうです。 背景をもう少しちゃんとフォローしてもう少しバリエーションのあるセットで試してみようと思います。 あとは配座発生も考えるとそれなりにバルキーな計算になりそうです。