align-itを使ってみた

以前からファーマコフォアベースで分子をアライメントするツールをささっとかけると
良いな〜って思ってました。
社内リソースを使おうと思うと、ライセンスの限りがあるのでオンデマンドに解析できなかったりするためです。
RDKitでやっている事例があったのでそれを参考に書いてみたのですが、途中でだんだん複雑になってしまい、ちょっとそのまま放置しています。
そこで、もう少し簡単にということでalign-itを使ってみました。
これは前にも紹介したsilicos-itから公開されているオープンソースです。
Openbabelがカバーするフォーマットを入出力に使えるので、特定のフォーマットに変換するという作業から解放されてよいと思います。
インストールはマニュアルに従って

> cd /usr/local/src
> sudo tar -xvf ~/Downloads/align-it-1.0.1.tar.gz
> cd align-it-1.0.1
> sudo mkdir build
> cd build
> sudo cmake ..
> sudo make
> sudo make install

でよいと思います。ubuntuにパッケージマネージャーでopenbabel入れてある状態でやると
openbabel2 無い。とエラーがでました。どうもパッケージマネージャーではな入らないみたいですのでこの場合は、自分でopenbabelを最初にビルドしてやるといいかと思います。

続いてテスト用のファイルの準備です。
今回は面倒なんでpubchemでEGFR阻害剤のデータ21化合物をとってきました。
二個目の分子をリファレンス(ref.sdf)としてアライメントします。

$ align-it -r ref.sdf -d pubmed_egfr.sdf -p out.phar -o out.sdf -s score.tab --rankBy TANIMOTO
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  Align-it v1.0.1 | Apr  1 2012 11:40:30

  -> GCC:         4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2336.1.00)
  -> Open Babel:  2.3.1

  Copyright 2012 by Silicos-it, a division of Imacosi bvba

  Align-it is free software: you can redistribute it and/or modify
  it under the terms of the GNU Lesser General Public License as published
  by the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  Align-it is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU Lesser General Public License for more details.

  Align-it is linked against OpenBabel version 2.
  OpenBabel is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation version 2 of the License.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

COMMAND_LINE OPTIONS:

  -> Reference file:    ref.sdf
  -> Reference type:    MOL
  -> Database file:     pubmed_egfr.sdf
  -> Database type:     MOL
  -> Mol output file:   out.sdf
  -> Pharm output file: out.phar
  -> Scores file:       score.tab
  -> Cutoff:            no
  -> Best hits:         no
  -> Rank by:           0
  -> Functional groups: AROM HDON HACC LIPO NEGC POSC HYBH HYBL
  -> Hybrids:           yes
  -> Epsilon:           0.5
  -> Merge pharm:       no
  -> Include normals:   yes
  -> With exclusion:    no
  -> Scores only:       no
  -> Quied mode:        no

Reference pharmacophore 5328042
   number of points:            8
   number of exclusion spheres: 0
   totalvolume:                 192.851
..
Processed 21 molecules
21 alignments in 0.849194 seconds (24.7293 alignments per second).

結果

$ more out.phar
NAME
HYBL    -1.24093        1.32284 0.407397        0.7     0       -0.120027       -0.867182       0.289608
HYBL    -3.03966        -0.157844       -0.170924       0.7     0       -0.120027       -0.867182       0.289608
HDON    0.564533        -0.760158       0.474673        1       1       0.236509        -1.68345        0.274832
HDON    -3.99885        -2.74223        -0.736352       1       1       -4.34047        -3.66037        -0.937119
HYBL    -1.01379        2.35455 0.601621        0.7     0       -0.120027       -0.867182       0.289608
HYBL    -3.67946        0.960047        -0.140495       0.7     0       -0.120027       -0.867182       0.289608
$$$$
NAME
HYBL    -1.23497        1.32782 0.409317        0.7     0       0       0       -1.11022e-16
HYBL    -3.03463        -0.151717       -0.169017       0.7     0       0       -1.66533e-16    -2.77556e-17
HYBL    3.35087 -0.639483       0.171317        0.7     0       0       0       0
HDON    0.5692  -0.7563 0.4762  1       1       0.208101        -1.67633        0.628354
HDON    -3.9954 -2.7354 -0.7347 1       1       -4.33756        -3.65332        -0.935583
HYBL    -1.04641        2.3288  0.59128 0.7     0       2.22045e-16     -4.44089e-16    -1.11022e-16
HYBL    -3.66258        0.968428        -0.135842       0.7     0       -4.44089e-16    0       2.77556e-17
HYBL    4.39451 -0.954693       0.404542        0.7     0       8.88178e-16     0       -5.55112e-17
$$$$

.pharにはフィーチャー ファーマコフォアの三次元座標 ガウシアン体積のα値 方向性があるかないかのブール値 そのベクトル
のタブ区切りがかえってきます。このデータって何で見れるんだろう???

またオプションの -sはスコアを返すので

lion:pubmed_egf iwatobipen $ cat score.tab
5328042	192.851	5328028	165.959	138.741	0	138.741	6	0.6304	0.7194	0.836
5328042	192.851	5328042	192.851	192.851	0	192.851	8	1	1	1
5328042	192.851	9818251	165.959	150.395	0	150.395	7	0.7216	0.7798	0.9062
5328042	192.851	9882519	192.851	163.02	0	163.02	7	0.7321	0.8453	0.8453
5328042	192.851	9885081	165.959	150.406	0	150.406	7	0.7217	0.7799	0.9063
5328042	192.851	10220590	219.743	191.497	0	191.497	8	0.8661	0.993	0.8715
5328042	192.851	10222656	192.851	178.099	0	178.099	8	0.8579	0.9235	0.9235
5328042	192.851	10245856	150.209	149.599	0	149.599	6	0.7733	0.7757	0.9959
5328042	192.851	10276061	192.851	163.018	0	163.018	7	0.7321	0.8453	0.8453
5328042	192.851	10868706	123.318	94.8864	0	94.8864	4	0.4288	0.492	0.7694
5328042	192.851	10870494	177.101	136.714	0	136.714	6	0.5862	0.7089	0.772
5328042	192.851	10902421	139.067	110.037	0	110.037	5	0.4959	0.5706	0.7912
5328042	192.851	10905749	165.959	158.397	0	158.397	7	0.7904	0.8213	0.9544
5328042	192.851	10926736	150.209	94.7072	0	94.7072	4	0.3813	0.4911	0.6305
5328042	192.851	10938977	192.851	150.2	0	150.2	7	0.6378	0.7788	0.7788
5328042	192.851	10992313	150.209	149.602	0	149.602	6	0.7733	0.7757	0.996
5328042	192.851	11012217	80.6759	79.1852	0	79.1852	3	0.4075	0.4106	0.9815
5328042	192.851	11055442	123.318	94.7115	0	94.7115	4	0.4277	0.4911	0.768
5328042	192.851	11080556	192.851	162.84	0	162.84	7	0.7307	0.8444	0.8444
5328042	192.851	11807783	165.959	163.157	0	163.157	7	0.8339	0.846	0.9831
5328042	192.851	11822927	123.318	109.048	0	109.048	5	0.5265	0.5654	0.8843

となりそれぞれの列はマニュアルにあるように
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 of the two structures
6 Overlap between pharmacophore and exclusion spheres in the reference
7 Corrected volume overlap between database pharmacophore and reference
8 Number of pharmacophore points in the processed pharmacophore
9 TANIMOTO score
10 TVERSKY_REF score
11 TVERSKY_DB score

となってます。
2列目と5列目が近い程よく重なってるということになります。
また -o 指定をすることでアライメント後の分子を書き出すことにしています。
CUIしか無いですが、結構有用なツールと思いました。

Advertisements

特許をパースする

不勉強にてxmlって
ブラウザで見られるけど何これ?
と以前思ってました。、、、恥ずかしい
よくよく調べるととても便利そうです。

で、xml形式の特許があると解析に使えそうと思うわけです。
集めてきた特許から、番号や、優先日、出願人とかの情報をささっと一覧にできると
いいよなーと思っていたこともあったので、

google先生に聞くとそんなことずいぶん前からなされているようでした。。。

さて、xml形式の特許が結構簡単に入手できる環境が整ったので
テスト的にスクリプトを作ってみました。
これは同じフォルダにある.xml形式のファイルをを全部読み込んで出願人とかをタブ区切りで出力します。
 フランス語とか、エンコでこけるのでアブスとは英語に限定しました。

相当行けてないスクリプト感ただよいますが、一応動いた。のでモチっと改良する予定

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from xml.etree.ElementTree import *
import sys, glob
patentlist = glob.glob("./*.xml")

f = open("pase_result.txt", "w")
f.write("PATENT_NO\tTITLE\tAPPLICANT\tABST\n")
    
for patent in patentlist:
    tree = parse(patent)
    root = tree.getroot()
    patent_no =  root[0][0][0][0].text+root[0][0][0][1].text+root[0][0][0][2].text
    title = root.findall(".//invention-title")[0].text
    applicants = root.findall(".//orgname")
    applicant = ""
    for  i in applicants:
        applicant += i.text+","
    abst = root.findall(".//abstract")
    for i in abst:
        if i.get("lang")  ==  "eng":
            l = i.getiterator("p")
            abstract = ""
            for i in l:
                abstract += i.text
            f.write("%s\t%s\t%s\t%s\n"%(patent_no, title, applicant, abstract))
        else: pass
   
f.close()

xmlの階層構造がそもそも勉強に使っていた例に比べるとかなり入り組んでいたのでまだ不具合たくさん出そうです。
がんばると、テキストなのでバラバラにして複数の出願人がどれくらいの頻度出るかとか、簡単に見れそうです。

最近読んだ文献

行き帰りの電車は文献を読むのにちょうどいい時間だったりします。昼はPythonのスクリプト書いたりRSS眺めてることが多いかなぁ。
 週末良さげな文献にほう見つけた。どっちも前に報告はあったけど今回のはそれがまとまってて自分的には刺激的。

A modern in vivo pharmacokinetic paradigm: combining snapshot, rapid and full PK approaches to optimize and expedite early drug discovery
 H to L , Lead optimizationの過程において、PKは重要度が高いが、なかなか試験ができないことがままあります(?)。
化合物の主活性ばっかり追い求めて、、、、 はて。と悩むのは行けてないわけですよねぇ。
 このレビューではsnap shot PK, rapid PK, full PKそれぞれについての比較を分かりやすくまとめてくれております。
 この前のsnap shot PKにフォーカスしたレビューでは全体のシステムの話が詳しく書かれており、こちらも勉強になりました。
 snap shot PK, rapid PKはSDがないので、個体差が大きい薬剤の場合多少ぶれるのかもしれません。
でもfull PKとだいたい良い相関示していると書いてあったので、そもそも化合物のクオリティーが高いのかもしれません。
 IT進んで情報処理系は結構技術向上してますのでデータをしっかり取る部分がボトルネックになる。現状に甘んじること無くそこをしっかり改善して行くというアプローチは、見習わねば。

A chemistry wiki to facilitate and enhance compound design in drug discovery
 あったらいいよね〜っていろんな人に言ってみても”?”という顔をされたことがここにあった訳で。
These web-based solutions are transparent,
always up-to-date, open for all users to view and contribute to,
interoperative and completely customisable.
誰もが見れて、すぐにアップデートされて、互換性があって、カスタマイズできる。ところにみそがあるように思う。
データーベースというと、何か結果を貯めとくみたいな意味合いが強くて、こういった仮説とか、フィロソフィーみたいなもの
をちゃんと残すシステムって発想持ってる人自分の周りにはあんまりいない。
 デザインを書き出して、in silico等を利用して優先順位をたてて、実際つくって、結果を解析。
まさにDMTAのサイクルにリンクしている。
 Fig5を見るとサイクルの加速具合が一気にあがっていることにこのシステムのすばらしさが伝わる。

Adoption outside the initial user-group was driven by word of
mouth and so the Compound Design Database found positive
response from many users.
ということで、この辺に何となく触発されたり
PythonベースのDjangoが最近使われていたりするそうで。

ちょっとこれはハイレベルすぎるが、何かデザインをまとめていけるようなシステムあるといいと思う。考えてみよう
いつかの会議のパワポ。とか、何かナンセンスじゃないすかね。

shape-itを使ってみた

2Dの分子の類似性を比較する場合、各種FP計算を実施してあれこれ考えるわけですが、
社内的都合にて3Dの場合、ささっと使えるツールが無い訳です。世の中的にはOpenEye社のROCSが使われることが多そうな気がします。文献出の引用例や、ユーザー会のプレゼンをみてもかなり良さげです。
お気に入りのRDKitにはShapeTanimotoメソッッドがありますが、アライメントは別途かける必要があるのでちょいと面倒でございました。
そんなおり、某ワークショップで、silicos-itのshape-itは?と教えていただきました(ありがとうございます)。
winではちゃちゃっとBuildするのは無理そうなんでまずはMacでやってみます。
shape-itのマニュアルを読んでみるとベースにOpenBabelを使っているのでOpenBabelがフォローするファイルフォーマットならなんでもOKという点も良い感じです。
Buildはマニュアルに従って,

> cd /usr/local/src
> sudo tar -xvf ~/Downloads/shape-it-1.0.1.tar.gz
> cd shape-it-1.0.1
> sudo mkdir build
> cd build
> sudo cmake ..
> sudo make
> sudo make install

で準備は完了です。
続いてテスト用のファイルの準備です。全部openbabelでやればいいという噂もありますが、こちらは自分のスキルの都合上RDKitで取りあえずのファイルを作ります。テストファイルはpubchemからもらったDYRK1Kのデータ。
2Dの分子から3D配座を発生させて取りあえずUFFで一応最適化しておきます。

from rdkit import Chem
from rdkit.Chem import AllChem
mols = [m for m in Chem.SDMolSupplier( "test.sdf" ) ]
mols3d = [  ]
for m in mols:
    m = Chem.AddHs(m)
    AllChem.EmbedMolecule( m )
    AllChem.UFFOptimizeMolecule( m )
    mols3d.append( m )

out = Chem.SDWriter( "mols3.sdf" )
ref = Chem.SDWriter ( "ref.sdf" )

for m in mols3d:
    out.write( m )
ref.write( mols3d[0] )
out.close()
ref.close()

全部の準備が終わったので
実行してみました。
ここでは先頭の分子をreferenceとして使います。

$ shape-it -r ref.sdf -d mols3.sdf -o out.txt -s score.txt
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  Shape-it v1.0.1 | Jul 20 2012 21:18:30

  -> GCC:       4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2336.1.00)
  -> OpenBabel: 2.2.99

  Copyright 2012 by Silicos-it, a division of Imacosi BVBA

  Shape-it is free software: you can redistribute it and/or modify
  it under the terms of the GNU Lesser General Public License as published
  by the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  Shape-it is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU Lesser General Public License for more details.

  Shape-it is linked against OpenBabel version 2.
  OpenBabel is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation version 2 of the License.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


COMMAND_LINE OPTIONS:

  -> Reference file:    ref.sdf
  -> Database file:     3dmols.sdf
  -> Output file:       out.txt
  -> Scores file:       score.txt
  -> Best hits:         no
  -> Scoring only:      no
  -> Extra iterations:  no
  -> Rank by:           Shape-it::Tanimoto
  -> Cutoff:            no
  -> Output reference   yes

..
Processed 29 molecules
29 molecules in 1.34573 seconds (21.5496 molecules per second)

結果

$ cat score.txt

dbName	refName	Shape-it::Tanimoto	Shape-it::Tversky_Ref	Shape-it::Tversky_Db	dbName	refName	Shape-it::Tanimoto	Shape-it::Tversky_Ref	Shape-it::Tversky_Db	overlap	refVolume	dbVolume
99380785	99380785	1	1	1	226.68	226.68	226.68
99380786	99380785	0.849	0.933	0.904	211.8	226.68	234.63
99380783	99380785	0.743	0.849	0.856	192.47	226.68	224.83
90945032	99380785	0.713	0.822	0.843	186.12	226.68	220.44
90945029	99380785	0.713	0.799	0.868	180.31	226.68	206.63
90945025	99380785	0.695	0.86	0.784	195.92	226.68	251.28
90945024	99380785	0.556	0.719	0.71	163.17	226.68	229.97
90945021	99380785	0.708	0.814	0.845	184.11	226.68	217.45
90945018	99380785	0.747	0.81	0.905	182.63	226.68	200.38
90945016	99380785	0.638	0.812	0.749	184.93	226.68	247.97
90945014	99380785	0.636	0.754	0.803	170.25	226.68	211.23
90945012	99380785	0.596	0.733	0.761	165.79	226.68	217.41
90944997	99380785	0.674	0.785	0.826	177.44	226.68	214.08
90944996	99380785	0.477	0.661	0.632	150.27	226.68	238.47
90944995	99380785	0.707	0.825	0.832	186.87	226.68	224.39
90944993	99380785	0.702	0.889	0.769	203.22	226.68	266.14
90944991	99380785	0.681	0.816	0.804	185.21	226.68	230.49
90944990	99380785	0.581	0.724	0.746	163.94	226.68	219.41
90944987	99380785	0.559	0.738	0.697	167.89	226.68	241.73
90944985	99380785	0.789	0.847	0.92	191.18	226.68	206.7
85239768	99380785	0.716	0.874	0.798	199.13	226.68	250.67
85239764	99380785	0.741	0.855	0.848	193.99	226.68	229
85239750	99380785	0.57	0.745	0.708	169.28	226.68	239.76
85239713	99380785	0.72	0.781	0.902	175.69	226.68	193.18
85239696	99380785	0.605	0.766	0.743	173.87	226.68	234.36
85239693	99380785	0.58	0.685	0.79	154.17	226.68	193.44
85239685	99380785	0.68	0.818	0.802	185.54	226.68	231.6
85239684	99380785	0.752	0.859	0.859	194.6	226.68	226.58
4239891	99380785	0.621	0.727	0.81	163.9	226.68	201.09

volumeだけの比較になりますが、ちゃんとアライメントもかけて出力してくれそうです。
背景をもう少しちゃんとフォローしてもう少しバリエーションのあるセットで試してみようと思います。
あとは配座発生も考えるとそれなりにバルキーな計算になりそうです。

サファリパーク

今日は家族でサファリパークに行ってきた。少し曇りがちで涼しく良い気候だった。

最初に歩いてぶらぶらと動物を見て、それから車でパーク内を見学。もっと暑さでぐてっとしてると思ったけどそれなりに動物が起きてて子供も結構楽しそうだった。

カピパラ

 

カンガルー

キッツカレー

最初あまり期待してなかったけど思った以上に自分も楽しめた。