04012013 chemoinformatics bioinformatics
僕のフィールドには一応Bioinformatics,Chemoinformaticsも含まれているので、文字列だけでどこまでいけるのかは非常に興味がある。まぁ物理法則を無視できなくなってくると文字列処理ではどうしようもないんだけどね。
3章のBWTまで読んだけど、3-4のBWTの性質と復元のところがさらっと流されていてさっぱり分からなかったので、検索したらわかりやすい説明を見つけた。
が、全体的には丁寧に解説されていて、分かりやすいと思うし、読んでて楽しい。
個人的にはSMILESみたいな構造情報を文字列にしたものをBWTで扱えないかなぁと考えている(chemoinformaticsに応用してみたい)。構造変換ルールにもイディオムみたいなのあるしね。MMPを高速に検索できても嬉しいだろうし、なんか使い道がありそうな気がするんだけど。
01122012 chemoinformatics Python
ローカルに構築したChEMBLにアクセスする場合には生のSQLじゃなくてSQLALchemyを使ったほうが楽だ。その場合にはMetaDataのrelflectをTrueで。
from sqlalchemy import create_engine, MetaData db = create_engine('mysql://user:password@localhost/chembl_14?charset=utf8&use_unicode=0') metadata = MetaData(bind=db, reflect=True) table = metadata.tables['target_dictionary'] stmt = table.select(table.c['organism'] == 'Homo sapiens') for r in stmt.execute(): print r.chembl_id, r.pref_name
ChEMBLって、生物種を正規化してないのね。パフォーマンスの問題かな？
Bioinformatics Programming Using Python01112012 chemoinformatics
システム構築してるはずが、ここのところ分析多めというか120%
ちょっと調べ物をしていてCanopy Clusteringという手法がmahoutに実装されていて面白そうだったので、論文読んでみた。
最初、自動で分割数選んでくれるのかーと勘違いしてたのだけど、荒くクラスタリングしてその後K-meansとかGreedy Agglomerative Clusteringに渡すだけらしい。分割数も荒いクラスタリングの中でAICとかBICを使って普通に求めるっぽい。
調べていく過程でIIR読めというような啓示があったが、積んであるのでそろそろ読もうかと思った。
NNがアツいらしい。活性予測なんかでも好成績をおさめているらしいのでちょっと興味がわいた。
Merck Molecular Activity Challengeはなかなかいい試みだなぁ。日本でもどこかがやらないかなぁ(リーディングカンパニーとか)
特に興味を惹かれたのは
いくつかの分野、特に人の認識に関係するような画像、音声、言語といった分野では既存の技術を大きく超えている感があります。
ドラッグライクネスというモヤモヤ感満載の指標がRo5から置き換わるのであれば素晴らしいことかもしれない。
24102012 chemoinformatics work
redmineはGitのようなVCSをセットにして使うのが当たり前なので、創薬系でもそのようなバージョン管理システムがあればいいなぁと思っている。
現状の化合物管理データベースは履歴管理をしないし、コミットログに対応するものも記録できないので、単なる構造データ管理としてしか機能しないが、化合物管理データベースのフロントにwikiかそれをモディファイしてコミットログとか履歴管理に対応する機能をもたせたものを用意すればやりたいことができそうだ。さらにAPIを用意してcliのラッパーとかpython,rubyバインディングを用意すればredmineと連携させやすいだろう。
logでログとかコミットツリーが出て
$ chemit log commit C0000123 Author: ChemistA <XXXXX> Date: Sun Aug 19 10:05:53 2012 +0900 improve solubility commit C0000122 Author: ChemistB <YYYYY> Date: Sun Aug 19 10:04:58 2012 +0900 reduce clogp fixes #24 commit C0000121 Author: ChemistA <XXXXX> Date: Sat Feb 4 09:02:58 2012 +0900 代謝安定性を改善 fixes #23
それから
高速文字列解析の世界――データ圧縮・全文検索・テキストマイニング (確率と情報の科学)
Mahoutイン・アクション
情報検索の基礎$ chemit diff C0000121 C000123
と打つと構造の差分とMCSが表示されようなものを作ってredmineと連携できれば素敵じゃないかなと思ったので、実証用のコードを書いてみようかな。
19102012 chemoinformatics
Hit Finding -> Lead Finding(H2L) -> Lead Optimization まで全て経験するとMMPは表層的というか現象論的だよなぁと思う。SBDDみたいな分子認識ベースを土台にしている人にとってはMMPは受け入れがたいほど荒いモデルだなぁと感じるだろうしねぇ。
MMPとFBDDの関係とかLEとMMP(activity cliff)の関係をよく考えてみると視点が開けると思うのだけど。
例えば
といったあたりを考えればいいと思う。
MMPが表層的であると理解できれば、外れ値（つまりcliff）はインシデントとして周知されるべきであり、インシデント駆動開発として回すことができるんじゃないかなと思っているし、分析する側としても自動的に予想外の出来事が収集されて楽しいだろう。
Redmineを創薬向けにカスタマイズすれば色々面白いことができると思うんだけど、興味ある人いないかな。
17102012 chemoinformatics work
創薬のためのメトリクスとかアジャイル開発手法の応用みたいな内容を話しますが、資料の中ではコンテクストや背景や動機なんかを端折ってるために、全体像がつかめないかもしれないので、参考になりそうなエントリをまとめておきます。
アジャイルとはなんなのか？
成果主義を例に出すまでもなく、組織を活性化するためには賃金による報奨は必ずしもいいとは限らない。承認欲求に訴えるためにはどうしたらいいか？
加えて暗黙知を形式知になりやすくするようなカルチャーをつくるにはどうしたら良いか？といったこともシステム設計として考えていく必要がある。
生産のアナロジーが使えそうな仕事または労働時間と生産性が比例するような仕事。参考書籍は沢山出ているので乱読すればよいでしょう。
ELNは単に紙の実験ノートを電子化すればいいってもんじゃなくて、もっとリポジトリへのアクセス端末としての設計を考えなきゃいけないんじゃないの？とか思う。 認証だとかPart11対応だとかはバックエンドのデータベースの問題であってそんなのイノベーションとは関係無いだろ？と思っているのでELN界とは距離を置いてる。
15092012 chemoinformatics q-chem
懇親会でpygamessじゃなくてもFacioで良かろうと言われたのだけど、対話環境よりも、スクリプトでやれたほうがいい場面が多いので書いておく。ちなみに僕はAvogadroユーザーです。
適当なsdfが手元になかったので
C CC CCC
というSMILESのリストを
babel -ismi carbons.smi -osdf carbons.sdf
とやってsdfにしたものを使いました。
sdfの分子群に対して計算したい場合、対話環境だとチマチマと作業しないといけないが、スクリプトだとforループを回せばよいし、結果がmoleculeオブジェクトで返ってくるので、そのあとのケモインフォマティクス的な取り回しが楽です。Gamessのアウトプット睨みながらコピペとか苦行だし、アウトプットの処理用のパーサー書くならinputから全てスクリプトで処理したほうがいいんじゃないかなと。
>>> import pybel >>> from pygamess import Gamess >>> g = Gamess() >>> mols = pybel.readfile('sdf', 'carbons.sdf') >>> for mol in mols: ... nmol = g.run(mol) ... print nmol, nmol.energy ... C -39.7265813363 CC -78.3052918313 CCC -116.885136991
ケミストのように化合物を個別のものとして一つ一つ丁寧に対応していくのか、情報論的にまとめて取り扱いたいのかとかそういうあたりで使い分ければいいのだと思う。オービタル見たい場合にはGUI必須だけど、エネルギーみたいな数字だけに関心があるのであればGUIは要らないかな。
28082012 chemoinformatics Python
論文読んでたら、コンベンショナルなMMPの求め方が載ってたので、似たようなやり方を考えてみた。
要するに網羅的に結合を切断していって
{Scaffold; [fragment, fragment, fragment]}
っていうハッシュを作っていく。fragmentは共通のScaffoldを持つので配列中の任意のフラグメントのペアがMMPということになる。
僕の切断ルール
CHEMBL327743の結果
('CC[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1CC(=O)*)/C=C/c1ccccc1)C', 'COc1cccc(c1*)OC', 'CHEMBL327743') ('CC[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1C*)/C=C/c1ccccc1)C', 'COc1cccc(c1C(=O)*)OC', 'CHEMBL327743') ('CC[C@@H]([C@@H](/C=C/c1ccccc1)N*)C', 'COc1cccc(c1C(=O)CN1CCCC[C@H]1C(=O)*)OC', 'CHEMBL327743') ('CC[C@@H]([C@H](NC(=O)*)/C=C/c1ccccc1)C', 'COc1cccc(c1C(=O)CN1CCCCC1*)OC', 'CHEMBL327743') ('CC[C@@H](C(NC(=O)[C@@H]1CCCCN1CC(=O)c1c(OC)cccc1OC)*)C', '*/C=C/c1ccccc1', 'CHEMBL327743') ('COc1cccc(c1C(=O)CN1CCCC[C@H]1C(=O)N[C@H](/C=C/c1ccccc1)*)OC', 'CC[C@H](C)*', 'CHEMBL327743') ('*c1ccccc1', '*/C=C\\[C@H]([C@H](CC)C)NC(=O)[C@@H]1CCCCN1CC(=O)c1c(OC)cccc1OC', 'CHEMBL327743') ('*OC', 'CC[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1CC(=O)c1c(*)cccc1OC)/C=C/c1ccccc1)C', 'CHEMBL327743') ('*OC', 'CC[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1CC(=O)c1c(*)cccc1OC)/C=C/c1ccccc1)C', 'CHEMBL327743') ('*CC', 'COc1cccc(c1C(=O)CN1CCCC[C@H]1C(=O)N[C@@H](C(C)*)/C=C/c1ccccc1)OC', 'CHEMBL327743') ('CC[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1CC(=O)c1c(OC)cccc1OC)/C=C/c1ccccc1)*', 'C*', 'CHEMBL327743') ('*C[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1CC(=O)c1c(OC)cccc1OC)/C=C/c1ccccc1)C', 'C*', 'CHEMBL327743')
どこで切断されたかどうかを明確にするために切断されたボンドにダミーアトム(atomicnumが0のアトム)をおいた。
canonicalSMILESで出力して.の文字でsplitした文字列をそのままキーに使っているが、openbabelのcanonicalizationアルゴリズムはMorgan Algorithmを改変したものを使っているのでおそらく大丈夫。
下のコードでは標準出力に書きだしたが、実際はデータベース化することになる。多分Mongo使う。
smi_string,id = clone.write(format="can")[:-1].split('\t') frag1, frag2 = smi_string.split('.') dbh.mmp.update({"core": frag1}, {"$push": {"frags": (frag2, id)}}, True) dbh.mmp.update({"core": frag2}, {"$push": {"frags": (frag1, id)}}, True)
こんな感じだろうか。
以下MMPを求めるのに使ったコード
import pybel import openbabel import sys input = sys.argv[1] def mkclone(mol): obc = openbabel.OBConversion() obc.SetInAndOutFormats("mol", "mol") molstring = obc.WriteString(mol) new_mol = openbabel.OBMol() obc.ReadString(new_mol, molstring) return new_mol def ring_member(rings, atom): for i, ring in enumerate(rings): if ring.IsMember(atom): return i return -1 def is_same_ring_member(rings, a1, a2): a1_idx = ring_member(rings, a1) a2_idx = ring_member(rings, a2) if a1_idx == -1 or a2_idx == -1: return False return a1_idx == a2_idx def detect_del_bonds(mol): rings = mol.OBMol.GetSSSR() del_bonds = [] for a in mol: if a.atomicnum == 6: for na in mol: n_a = na.OBAtom b = a.OBAtom.GetBond(n_a) if b != None and b.GetBondOrder() == 1 and not is_same_ring_member(rings, a.OBAtom, n_a): if a.OBAtom.GetIndex() < n_a.GetIndex(): del_bonds.append((a.OBAtom.GetIndex(), n_a.GetIndex())) return del_bonds for mol in pybel.readfile("sdf", input): chembl_id = mol.data["CHEMBL ID"] print chembl_id del_bonds = detect_del_bonds(mol) for a1, a2 in del_bonds: clone = pybel.Molecule(mkclone(mol.OBMol)) # !!! index starts 0 but GetAtom starts 1 !!! ob_a1 = clone.OBMol.GetAtom(a1 + 1) ob_a2 = clone.OBMol.GetAtom(a2 + 1) bond = ob_a1.GetBond(ob_a2) clone.OBMol.DeleteBond(bond) new_atom1 = clone.OBMol.NewAtom() new_atom1.SetAtomicNum(0) new_bond1 = clone.OBMol.NewBond() new_bond1.SetBegin(ob_a1) new_bond1.SetEnd(new_atom1) new_bond1.SetBondOrder(1) new_atom2 = clone.OBMol.NewAtom() new_atom2.SetAtomicNum(0) new_bond2 = clone.OBMol.NewBond() new_bond2.SetBegin(ob_a2) new_bond2.SetEnd(new_atom2) new_bond2.SetBondOrder(1) smi_string, id = clone.write(format="can")[:-1].split('\t') frag1, frag2 = smi_string.split('.') #print (frag1, chembl_id) #print (frag2, chembl_id) print (frag1, frag2, chembl_id)
23082012 chemoinformatics
そのうちMMPベースのプロジェクト診断ツールでも作ろうかなと思っているので、fmcsを使ったMMP構築をやってみた。
fmcs sample.sdf --output-format complete-sdf --save-smiles-tag \ MCS_SMILES --save-atom-indices-tag MCS_ATOMS --output-all
で、MCSのSMILESと対応するatom indexをタグ付けしたsdfを出力しておいてpybelでMCSを削除する。
import pybel import sys input = sys.argv[1] for mol in pybel.readfile("sdf", input): print mol.data["MCS_SMILES"] del_atom_idxs = [int(num) for num in mol.data["MCS_ATOMS"].split()] del_atoms = [] for i,atom in enumerate(mol,0): if i in del_atom_idxs: del_atoms.append(atom) for a in del_atoms: mol.OBMol.DeleteAtom(a.OBAtom) print mol
そうすると、MMPが作れる。
もう少し具体的な例をと思ってChEMBLからsdfでも取ってきてとか思ったけど、やり方がわからなかったのでわからんとか言ってたら@32nmに教えてもらった。ありがとうございます。
というわけで、ChEMBLからカテプシンK阻害剤の化合物データからCHEMBL100616,CHEMBL101160,CHEMBL316685,CHEMBL327743,CHEMBL329731,CHEMBL419891の6化合物を引っ張ってきてやってみた
CCC(C)C(C=Cc1ccccc1)NC(=O)C1CCCCN1CC(=O)c1ccccc1 C(=O)(N)C CCC(C)C(C=Cc1ccccc1)NC(=O)C1CCCCN1CC(=O)c1ccccc1 Cl CCC(C)C(C=Cc1ccccc1)NC(=O)C1CCCCN1CC(=O)c1ccccc1 [NH+](=O)[O-] CCC(C)C(C=Cc1ccccc1)NC(=O)C1CCCCN1CC(=O)c1ccccc1 Cl CCC(C)C(C=Cc1ccccc1)NC(=O)C1CCCCN1CC(=O)c1ccccc1 [NH+](=O)[O-] CCC(C)C(C=Cc1ccccc1)NC(=O)C1CCCCN1CC(=O)c1ccccc1 OC.OC
あとはアッセイ情報と付きあわせてpIC50みたいなプロパティの引き算をすればMMPを利用した解析ができる。
ある程度コントロールできてる化合物群ではうまくいくんだが、MCSだと期待している以上にスキャフォールドをアサインしまうために、結果としてフラグメントが細切れになったり(Rの環が分割されてしまったり)してMMPを構成するにはちょっとあれかな。
21082012 chemoinformatics Python
openbabelはMCSを求められないのが難点なんだが、書けばいいよなと思いつつ実装してない。そうこうしているうちに、fmcsを見つけて、これ使えばいいかなーと思い始めた2者だけじゃなくて複数の化合物のMCSを求められるのも自分の求めていた機能だし。
fmcsを動かすにはRDKitをインストールする必要があって、「RDKitのインストールめんどいなー」とぶつくさ言ってたら@iwatobipenにsilicosを参考にしなはれと教えてもらった。
brewで入れられんのね。
> brew tap edc/homebrew-rdkit > brew install rdkit
これでサクッと入った。
fmcsのほうは
Redmineによるタスクマネジメント実践技法hg clone https://bitbucket.org/dalke/fmcs
cd fmcs
sudo python setup.py install
で入る。
さて、MLによると、
This work has been funded by Roche, with the explicit hope that it be incorporated into RDKit.
ということなのだが、確かRSOAPはPfizerが出してたはずだし、海外はopensourceとかそういうのに理解が深いんだろうか？また、そういうのに積極的な日本の製薬企業はあるのだろうか？なんて思ったりした(自分の観測範囲ではゼロなんだが)。
そういうことも含めて貢献できるような社会がくればいいのになぁと思う。
Design and Use of Relational Databases in Chemistry