久々にケモインフォクックブックを更新した。

先日同じことをするのにrxnフォーマットを扱えるライブラリ使ってちゃちゃっとやったけど、openbabelだとどう書くのがいいのかなと。

さて、SBDDとかLBDDとかそういう仕事に関わっているヒトはプログラミングくらい出来そうなイメージがあるけど、全然そんなことなくて、ちょっとしたシェルスクリプトも書けないヒトが多いです。基本アプリケーションの操作という作業がメインなので、プログラミング出来なくてもなんとかなるという。あとPipelinePilotみたいな便利ツールもあるしね。

という状況なので、仕事に関わりながらプログラミングを覚えていくというキャリアパスのようなもの存在しないと考えてよいです。それは単純に独学と自己啓発のエリアですな。

しかし、アプリケーションのみでは出来ない作業っていうのは存在するわけで、そうなるとやっぱプログラミングしなきゃいかんわけですよ。そこで、プログラミングが不得手なヒトは仕事の効率がガタ落ちするわ、最悪成果が出せないとか困ることになる。

僕なんかが「誰でもできるツールで誰でも出来る仕事をしてもしょうもないでしょ？」とかいうと「そういうツールを組み合わせて独創的な成果を」と反論されるんだけど、そんなこといってもしょうがないような。

だってそれは誰でもできるツールを使うヒトがごく自然にたどり着く結論というかみんなその制限の中で工夫しているからね。

そういった制限から抜け出すにはプログラミングは一つの解なのだろうなと思う。

Design and Use of Relational Databases in Chemistry
T. J. O'donnell
Crc Pr I Llc / ￥ 11,398 ()
通常2～5週間以内に発送

ま、とはいうもののプログラミングスキルは基本的に評価されないよね。

この本が欲しいが、ちょっと手が出せない

Optimizing the "Drug-Like" Properties of Leads in Drug Discovery (Biotechnology: Pharmaceutical Aspects)

Springer / ￥ 20,252 ()
通常2～4週間以内に発送

リードオプティマイゼーションってのは割と抽象的な言葉だなぁと思っていて、中期的にはターゲットタンパク質のオプティマイズ+physicochemical propertyの調整ではなくて、PharmacokineticsやPharmacodynamicsの最適化を指す言葉になるんだろうなぁと思う。

創薬ターゲットとしてのタンパク質は有限だって示されたわけなので、ある程度進んだあとの薬剤としての差別化ポイントは効き目の強さというよりは、効果の持続性とか投与経路が患者にとってより良いかとかのQOLでしょうしね。

そうすると臨床系のデータを探索の戦略に組み込めないとあかんよねとか。

javascriptで書かれたケモインフォ用のツールキット。

レンダラももちろんついてるが、firefoxだと表示されなかったのでSafariで。

Wikiみたいに適当な記法を用意してsmiles貼りつければ絵描いてくれるような仕組みにもできるので、これは色々と使い道がありそう。

• Cheminformatics in JavaScript: JChemHub
• jchemhub

最近ダラけたというか、やる気が下がっていたところ、stochastic proximity embeddingとつぶやかれていたので、おー面白そうと調べたらpdfがあったので読んでみた。

多次元尺度法同様に、距離情報から座標を構成する手法らしい。ペーパーだと次元縮約の方法として紹介されているけど、一回距離情報を求めてそれを任意の次元に置き直すのでまぁ似たようなものかなと。

ただ、SPEのほうはデータセットが大きくなっても計算量が爆発しないので大きいデータに使えるそうだ。

• マップする次元を決めたら初期値としてランダムに座標を与える
• ランダムに二点を選ぶ(xi, xj)
• 距離が一致するように遠ければお互いを近づけ、短ければお互いを離す

• 何回か二点を選んで更新処理を行ったらラムダの値を小さくして、上を又繰り返す

ほーこれで上手くいくんかいなと思ったのだが、プラクティカルにはよさそうかも。でも、ベストな構造におちないのと、ラムダの刻み幅は小さくしないといけないっぽいな。 ちなみに、CRANにもspeというパッケージがあったので、今回これを動かしてみた。

library(spe)
data(phone)
embed <- spe(phone, edim=2, evalstress=TRUE)
plot(embed$x)

距離行列から計算できると便利だけどこのパッケージは座標をインプットにしなきゃいけないのとユークリッド距離固定な感じですね。実用で使う場合は自分で実装した方が思い通りに動かせるかな。あとアルゴリズムが単純なので並行処理もできるような気がするけど。そのうちScalaで書いてみたい感じ。

ドッキングシミュレーション

というシュールなネタを思いついた。

参考

• paper child

朝から、openbabelのjavaインターフェースを使えるようにしようといじっていたのだけど、java.lang.UnsatisfiedLinkErrorを解決することができないまま夕方になってしまった。

んーCDK使うことになるのかなぁと思ったがMXのほうが楽そうなので今回はこっちを使った。

import com.metamolecular.mx.calc.MassCalculator
import com.metamolecular.mx.io.daylight.SMILESReader

object MxTest {
  def main(args: Array[String]) {
    val calc = new MassCalculator
    for (mol <- args) {
      val mw = calc.findAveragedMass(SMILESReader.read(mol))
      println(mol + ": " + mw)
    }
  }
}

実行してみる

scalac mxtest.scala 
scala MxTest "CCC" "CO" "CNC"
CCC: 44.09562
CO: 32.04186
CNC: 45.08368

Liftを使って構造関係を扱えるWikiを作ってみたい。

Scalaスケーラブルプログラミング[コンセプト&コーディング] (Programming in Scala)
Martin Odersky,Lex Spoon、Bill Venners
インプレスジャパン / ￥ 4,830 ()
在庫あり。

ほとんど化合物情報をTokyo Cabinetで管理してみると同じノリで出し入れできそう。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

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

import pybel
import couchdb

sdffile = "pc_sample.sdf"
couch = couchdb.Server()
db = couch.create('pubchem')

mols = pybel.readfile("mol",sdffile)

for mol in mols:
    title = mol.OBMol.GetTitle()
    molstring = mol.write("mol")
    db[title] = {"mol":molstring}

化合物情報みたいな、RDBで管理しやすいようなデータよりは、in vivoの薬理試験とか動態試験みたいな、プロトコル間で比較があまりなくて、かつ所見とかの文章データが重要なもののほうが向いているのかもしれんなぁと思った。PKデータなんてIDで探して、時点と個体のデータがJSONで引っ張ってこれるようにしとけばナイスすぎる！

javascriptからRの関数呼べんかな？そうすればJSONでデータ受け取ってPKfitとか使えんのにな。それかbrewとか

書籍もいくつか出てるけど、日本語ないなぁ。そのうち出るんかな？

CouchDB: The Definitive Guide
J. Chris Anderson
Oreilly & Associates Inc / 3390円 ( 2010-01-15 )

Beginning CouchDB
Joe Lennon
Apress / 3783円 ( 2009-12-30 )

Couchdb in Action
Christopher Chandler
Manning Pubns Co / 4180円 ( 2010-12-28 )

• CouchDB: The Definitive Guide
• CouchDB: the Definitive Guide 日本語訳

ふと、つぶやいた。

例えば

• プロジェクトで500検体合成
• ひとつの親あたり5-25化合物くらいの子を合成する
• プロジェクトの人員の制限により、一度に10ラインまでしか同時に走らない
  • 11番目以降はそれ以降のオプティマイゼーションは行われない

その時知りたいことは、

• {計画、合成、アッセイ、解析}というサイクルは何回まわるのか、
• 有望そうな化合物だけど、人的リソースの関係で選択されなかった歴史はどこにあるか？
• バックアップ候補の化合物はどこら辺で現れるか（初期？中期？）

あたり。

現実の系に近づけるために、親化合物にはポテンシャルがあって、近傍探索すると、ポテンシャルの近くで活性が変化するモデル。あと同じライン（ブランチ）を継続して合成すると、ポテンシャル曲線が底に近づくのと、合成もSARに関する知識がついてくるので、活性は上がりやすくなって、かつ変動幅も小さくなるようにweightはだんだん小さくなるようにしてみた。

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

# kzfm <kerolinq@gmail.com>

from random import random
import networkx as nx
import matplotlib.pyplot as plt


class Branch(object):
    """LeadOptimization flow"""

    count = 0

    @classmethod
    def inc_count(cls):
        cls.count = cls.count + 1
        return cls.count

    @classmethod
    def get_count(cls): return cls.count

    def __cmp__(self, other):
        return cmp(self.act, other.act)

    def __init__(self,potency,weight):
        self.num = Branch.inc_count()
        self.potency = potency
        self.weight  = weight
        self.act = self.potency + self.weight * random()

    def make_child(self,num_childs,potency,weight):
        return [Branch(self.potency + self.weight * random(),self.weight * 0.7) for i in range(num_childs)]

if __name__ == "__main__":
    max_comps        = 500
    initial_branches = 3
    nodesize         = []
    nrange           = []
    erange           = []
    generation       = 1
    heads = [Branch(0.3,1) for i in range(initial_branches)]
    G=nx.Graph()

    for h in heads: 
        nodesize.append(h.act * 30)
        nrange.append(generation)

    while Branch.get_count() < max_comps:
        new_heads = []
        generation += 1
        for branch in heads[:10]:
            for new_comp in branch.make_child(int(5+20*random()),branch.potency,branch.weight):
                nodesize.append(new_comp.act * 30)
                nrange.append(generation)
                erange.append(generation)
                G.add_edge(branch.num,new_comp.num)
                if new_comp.act > 0.75:
                    new_heads.append(new_comp)
        heads = new_heads
        heads.sort()

    nx.draw_spring(G, node_size=nodesize, node_color=nrange, edge_color=erange,alpha=0.7, \
cmap=plt.cm.hot, edge_cmap=plt.cm.hot, with_labels=False)
    plt.savefig("proj.png")

んーなんかいまいち。一応活性の強弱にあわせてノードのサイズが変わるようにしてみたんだけど、大きさがわからん。

Pharmacokinetic-Pharmacodynamic Modeling of Biomarker Response and Tumor Growth Inhibition to an Orally Available cMet Kinase Inhibitor in Human Tumor Xenograft Mouse Models という論文が良かった。というか、最近ボトムのあたりを這っていた僕のモチベーションがかなり上がった。

ケモインフォマティストとして、このレベルのPK-PD解析を出来るようにしたい。

この本の23,24章に最適化フェーズにおけるPK-PDモデリングとかPBPKモデリングの意義とかがあって、その内容が良い。

Preclinical Development Handbook: ADME and Biopharmaceutical Properties (Pharmaceutical Development Series)

Wiley-Interscience / ￥ 18,701 ()
通常2～4週間以内に発送

そして、ケモインフォマティストなのでSASとかNONMEMではなくてRを使う。

「非線形混合モデルを理解するために」というpdfも見つけた。

単にIDをキーにしてsdf呼び出すだけだったらKVSで良くね？と思ったところ、そういえばTokyoCabinetのPythonバインディングがあることを思い出したので、ちょっと試してみた。

データはこの時に落としておいたsdf

import pybel
import pytc

sdffile = "pc_sample.sdf"
db = pytc.HDB('compounds.db', pytc.HDBOWRITER | pytc.HDBOCREAT)
mols = pybel.readfile("mol",sdffile)

for mol in mols:
    title = mol.OBMol.GetTitle() 
    db[title] = mol.write("mol")

100万化合物くらいなら全部メモリに載せられそうな感じ。そしたら速いよねー、CouchDBも気になってきた。

Software Design ( ソフトウェアデザイン ) 2010年 02月号 [雑誌]

技術評論社 / ￥ 970 ()