blog of personal interest
RDKitにある２Dファーマコフォアを使って シミラリティマトリックスを作るスクリプトを書いてみた。 下のデフォルト設定では3000ビット以上のFPになるので 適当なbinに指定し直した方が計算が速くなるかも。 分子によってはGitBitIdxでエラーになるのでそれを回避するため try/except処理にした。エラーが出る分子については後で原因を追う。 TCの計算はお決まりの書き方でやると、エラーにうまく対処できない。 見にくいけどdict型でindxとFPの関係を保つようにしてみた。 よくよく考えると、SVMとかするならFPをbitstringで返しておしまいにした方がいいと気がついた。 使い方 $ python p4core.py で入力SDFに関してシミラリティマトリックスを出力される。 ECFP/FCFP/MACCSのTC 諸々と合わせて出力させると違いが分かるかな。
今日は昼休み中に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表記が未だにしっくりこないのでその辺が課題です。 ちょっとコード自身もきれいじゃないですね。
PubMedやCHEMBLはRESTFULサービスを提供しています。 その前は多分SOAPサービスでした。 Pubmedで検索かければいいと言われそうだけど、何か面白いことできないかと 思っておりましたら、ちょっと惹かれるものを発見。 PubChemPy 早速使ってみましょう 例に従っていきます csはリストオブジェクトなので 要素を取り出して とかできます。 と３次元のクエリにすればエネルギーもとって来れます。この辺はPubChem3dに詳しく記載があったはず。 aidからアッセイ情報を取得する場合はget_assayでとれるので とリスト型をうけて にてアッセイ情報をとりにいけます。 アッセイ情報にターゲットクラスも入っているのでヨサゲ。 データ自体はSWISSPROTからきているのかな u’Db Source: SWISS-PROT’, u’Db Version: 2012_01′, u’Multi: No’, u’Complex: No’, u’Relationship Type: Direct protein target assigned’, u’Target Class L1: Enzyme’, u’Target Class L2: Kinase’, u’Target Class L3: Protein Kinase’, u’Target Class L4: Tyr’, u’Target Class L5: Tk’, u’Target ClassContinue reading “PubChem rest”
pythonからインターネットにアクセスするパッケージはurllib2,requestsなどをつかいますが プロキシ経由でアクセスする必要がある場合も考えねばなりませぬ。 認証が無いときは というようにopenerを設定して開けますが、ベーシック認証が挟まると これではつながらないので とauthハンドラを入れておけばいい。 だけどコードの引数に直接ユーザー名とパスワードを入れるのが気持ち悪いなあと感じていたら、getpassというライブラリがあった。 これはプロンプト上にpass入力を求めてかつ、入力は隠される。 これだ！ で、 getuserメソッドはログインユーザーのユーザー名をとってくるのでAD環境なら これでIDは引っ張れる。 getpass内のプロンプトは適当。 これでいい感じになった。 けどまだまだスキル不足 もっといい方法があれば是非教えてください。
最近Rの勉強をしています。 TSは散布図は一度に一個なので、フィールドが複数ある場合不便です。 その点Rはデーターフレームであればplotで一気に描けるのでいいなあって思います。 そんなRできれいな図を描くツールとしてggplot2があります。 さらにプロットマトリックスを作るツールとしてGGallyが素敵です。 インストールはこんな感じで でggpairsでマトリックスプロットを作ります。 data引数にデーターフレームを columnsに使用するカラムを colors にファクターを入れると とすると、以下のようなグラフができます。 いいでね
RDKitとPythonを使ってデータを加工して Rで解析するというケースを頻繁にやっていると全部つながるといいなーと思います。 RPY2を使ってうまく繋げるのもいいですがKNIMEを使うとそれぞれの連結ができて便利かも。 ということでトライしてみました。 KNIME上のRノードでは全部Rという名前のオブジェクトに対して処理をするようで それを理解するのに時間がかかりました。 例えばPLS回帰をする場合はこんな感じです。 Rノードおのおのでplsパッケージ読み込まないと、 おいおいmvrなんて扱わないぜ。とおこられます。 で作ったのはこんな感じです。 一個はモデルを作ってモデルと実測を比較するため もう一個はモデルに当てるためのフローです。 RDKITノードでディスクリプタは簡単にとれるし、そのあとのモデリングは上のようなスニペット書いとけば後はよろしくな感じでいいので、楽できそうです。 ローカルで全部やってるからパフォーマンスが悪い、、、 まあいいか。
IUPAC, CAS_Nameから構造を起こすツールは手元にあるのですが、 ケミストとの視点で考えると使いにくかったり、一般名から構造が起こせると言いな。というニーズはあると思います。 この辺りはOPSINを使うのがいいと思われます。 がとあるところからCIRpyというクールなツールを見つけました。 Chemical Identifier Resolver(CIR)のPython インターフェイスだぜ！ということです。 OPSIN, ChemSpider当たりも入っているぜぃ。 ということで何か惹かれます。 ソースがGitHubにあがっているのでこれをゲットします。 使い方は直感的に行ける感じです。 という感じです。 CASナンバーもさっくりと。 その他に stdinchi stdinchikey inchi smiles ficts ficus uuuuu hashisy sdf names iupac_name cas chemspider_id mw # Molecular weight formula と色々なアウトプットが行けます。 それ以外に、記述子も行けます。 h_bond_donor_count h_bond_acceptor_count h_bond_center_count rule_of_5_violation_count rotor_count effective_rotor_count ring_count ringsys_count 構造式を落とすこともできますね。 これでカレントディレクトリにsdfが落ちます。素敵ですね。
web.pyでファイルをユーザーから受け取ってxlrdでDB登録用のフォーマットに直したい。 と思ったんですが、これがなかなかうまく行かないので難儀しました 結局open_workbook()の引数にweb.input()をそのまま渡す場合は file_contens = web.input(myfile={}).read() で行けることが分かりました。 最初は閉じたfoutをとりにいったのですが全然うまく行かない、、、参りました。 やりたいのはここで吸い出したデータをMS Accessに突っ込むことだから 今度はpyodbcを見るかなあ。 でも社内環境だから渋々の選択だが もっと別のDBの方がいいな。
web.pyを使ってサービスを提供できるようになったので 何か使えないかな〜と考えて見た。 任意のSDFがあったらPCAやって結果を返すとかあると PCAってなに？みたいな人もspotfireで図が見えるしいいかもと思った。 PCRの方がより良いと教えていただいたのでこちらの実装も考えよう。 さて、自分がやるならRを使うのだが、web.pyに入れるとなると全部pythonがいい。 rpy2入れたしこれで。と思っていたら動かないで肩すかし、、どうも社内の環境のwinで構築したR2.10はだめなようだ RDKitのモジュールはもうメンテしてないし、matplotlibがいいじゃないと言われたのでそれで取りあえず書いてみました。 化合物数＞ビット（2048）の場合はそのままでいいみたいですが、 化合物数＜ビットの場合はそのままではエラーになるので転置します。 で最後にもう一回スコアを転置して第三成分までを出力。 手ものとのデータでやったらあまりぱっとしない感じ。 もちっとデータ集めてやってみる。
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も理解せねば こういったこと、社内で言ってもきっと理解されないなぁ 何か考えるかな〜
社内はなんだかんだでエクセラーが多いので 基本エクセルベースのデータがたまります。 解析しようにもデータが正規化してないのでまずその加工から入るのですが ルーチンだしモチベーション上がりません。 マクロ書くにもVBA勉強すんのもなーと思っておりましたら、pythonで行けるという事を教えていただきました。（感謝） ググるとどうもxlrd, xlwd, xlutilsが良いようです。xlsxもサポートしているので社内の最近の環境でもバッチコイ。 でこれらはeasy_installかpipでさっくり入ります。winならインストーラーもあるようです。 早速使ってみました。 とりあえずやりたいのは、 ID ARGET_A TARGET_B TARGET_C TARGET_D cmp1 1000nM 100nM 1000nM 500nM cmp2 600nM 100nM 1050nM 500nM みたいなやつを ID TARGET MARK VALUE UNIT cmp1 A = 1000 nM cmp1 B = 100 nM cmp1 C = 1000 nM cmp1 D = 500 nM みたいにしたいのです。 マクロ書けば簡単でしょ。とかいう人いますけど、言ってる人が書かない。 で、ターゲットとあるA、B、CはCYPでもHL60でも薬理活性でもいいですね。Continue reading “pythonからエクセルのデータを加工する”
WIPOのpatentscopeで作成した RSSからタイトルとサマリを抜き取って 英数字以外で区切った単語リストを作るスクリプト。 sys.argv[1]に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を書いて 取りあえず動く。 リダイレクションをファイルにして ということで出力先をもう少し整形しようと思う。 プロキシハンドラも入れないとだめですね。
文献や特許関連の情報はRSS登録していると便利ですね。 ACS系のジャーナルはグラフィカルに見えるので好きです。 集合知プログラミングという本にfeedparserというライブラリの使用例がのっていて 便利そうでしたので使ってみました。 対象はjmc lettです。 summaryで全然情報がとって来れない。 アブストの情報をとるのはもう少し工夫が必要そう。 この辺のことをやるのはwebの知識も必要なのかな
各分子の最安定構造のモーメントと 発生させた配座郡のシミラリティ計算になるようにした。
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のデータで何か解析してみないと。
SDFには最小のエネルギーを持つ配座だけを テキストには全部のログを残すようにしてみた。 ちゃんとテキストの方には100分子で１０配座のデータがきました。 時間がちょっとかかるのでもっと荒くていいならiter =200くらいでいいかなぁ。
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)で行けそうだから後で追加しよう。
職場でRを使おうとするとライブラリのインストールに四苦八苦する。 プロキシを超えねばならないのです。 調べると、Rのショートカットの後ろに –internet2と入れたらいいとあるのですが これだと、ブラウザのプロキシ設定は読むけど、認証は別途入れる形なのでリポジトリにつながらない。 さらに調べた結果、 ホームディレクトリに .Rprofileを作って Sys.setenv(http_proxy=”http://proxy:port/” Sys.setenv(http_proxy_user=”ask” と入れると良いと分かった。 ちなみにこれを書いていても–internet2の記載を残しっぱなしだとそっちの設定になるみたいで つながらない。 このことに気がつくのにかなりの時間を費やした。。。。 今はいい感じになったのでまずはe1071, rcdkを入れました。
以前からファーマコフォアベースで分子をアライメントするツールをささっとかけると 良いな〜って思ってました。 社内リソースを使おうと思うと、ライセンスの限りがあるのでオンデマンドに解析できなかったりするためです。 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を使ってみた”
Is life worth living? 