献本ありがとうございます

Rによるバイオインフォマティクスデータ解析　第2版 －Bioconductorを用いたゲノムスケールのデータマイニング－
樋口 千洋
共立出版 / 4830円 ( 2011-01-08 )

内容はバイオインフォマティクスに限らずに割と広い内容をカバーした感じで、クックブックと逆引きの中間的なスタイルと言えば良いのだろうか？

Rのインストールから基本的な操作は（大体どの本にもあるように）載っていて

データマイニングとしては

• PCA
• ICA
• PLS
• MDS
• SPE
• k-means,Fuzzy cmeans
• spectral clustering
• NMF
• SOM
• decision tree
• kNN
• SVM
• RF
• LASSO
• MARS

がサンプルコードとともに簡潔に説明されている。

8章はバイオ系データの解析、チップとか。odesolveを利用したシミュレーションのサンプルもあって、SBMLRは面白そうだなぁと思った。メカニズムがどうなっているのかはモデルと実験系の不一致をよく突き詰めて考えることでしかきちんとした理解は得られないと思っている。

最後のほうの章は統合環境、データベースの連携、サーバー構築あたりの話。

あと、twitteRを使って生体組織名をつぶやくとその組織の図が返ってくるという例が載ってた。ちなみにRでもPitがあります。

ボロノイ図おもしろい。IBIS行ったときに、最近傍探索はボロノイ図を使えば効率的に求めることができると聞いて興味をもった(Blogopolisから学ぶ計算幾何にも載ってた)

本書は、実装が載っていないのが残念だけど、実例を交えて進んでいくので非常に分かりやすかった。

なわばりの数理モデル -ボロノイ図からの数理工学入門-
杉原 厚吉
共立出版 / ￥ 2,730 ()
通常2～5週間以内に発送

• 局所ドロネー性
• ドロネー辺は与えられた点集合に隣同士という関係を定義するものだと解釈できる
• 最大空円
• ロイド法
• 最遠点ボロノイ図と最小包含円
• 最小全域木の辺は相対近傍グラフの辺

ケミカルスペース(PCA)に対してボロノイ図描かせてPCRすれば、ケミストにとって分かりやすい(既存のどの化合物に近いか一目でわかる)表現方法にならないだろうか？

これは購入しておく。

Free Energy Calculations in Rational Drug Design

Springer US / ￥ 13,115 ()
一時的に在庫切れですが、商品が入荷次第配送します。配送予定日がわかり次第Eメールにてお知らせします。商品の代金は発送時に請求いたします。

Free energy calculations represent the most accurate computational method available for predicting enzyme inhibitor binding affinities. Advances in computer power in the 1990s enabled the practical application of these calculations in rationale drug design. This book represents the first comprehensive review of this growing area of research and covers the basic theory underlying the method, numerous state of the art strategies designed to improve throughput and dozen examples wherein free energy calculations were used to design and evaluate potential drug candidates.

うちの職場は古きよきというか（略）

量子化学計算のコストも下がったというか、小さい化合物だったら普通に待てる時間でレスポンス返せるのでWebインターフェースつくってある。配座の安定性調べたり、ソルベーションエナジー計算したりとか色々やれるので、合成力(how to make)よりも、何を合成するか(what to make)？つまり、どういう仮定を置くか、論理をどう組み立てるか？という化学反応の背景にある物理化学とか量子化学に対する理解を強めるのは重要でないのかなぁと思っているんだけど。

どうかな？

量子化学計算できないケミストは、BLAST検索できない自称分子生物学者みたいなもんじゃねーの？

と常々思っている。

で、そういうヒトってなにに基づいてケミストリーしてんだろ？という素朴なギモンがあるのだが、ここではそういう部分は抜きにして、とりあえずBLASTのようにweb上から簡単に計算できるような量子化学計算サービスを用意すれば、リテラシーの問題で量子化学計算に手を出せなかった層に訴求できるんじゃないかな？といつも思っていたのだけど、所詮他人事なのでスルーしてきた。

でも、いろいろスルーできない状況になってきたので、余暇を利用して少しそこら辺も考えてみたというか、ハッカソンのテーマとして自分に課したら、結構はかどったので良かった、素晴らしい。

というわけで、ケモインフォクックブックを更新した。

次回はIzu.R #1として、発表付きのハッカソン形式としてやれればいいかなと思っている。

熱海でだらだらと温泉につかりながら、コード書いたり書かなかったりという場を設けました。キーワードはchemoinformatics,bioinformatics,R,Python,Rubyで製薬企業でコード書いているヒト多めです。僕はFlaskいじるかopenbabelのMCS実装かなんかをする予定です(でもいい日本酒がゲットできたら飲んでるかも)

• 2010.10.29-30
• 芳泉閣@熱海
• 一泊二日で10500くらい

まだ何人か分の空きがありますので興味があれば私にメールかtwitterでメッセージを下さい。金曜から土曜にかけてという日程です。

今のとこの参加者

• kzfm
• zgmfx20a
• bohohu
• garuby
• shirahakase
• hiro_h

久々の更新だ。

• アロマティックなヘテロ環を漏れなくさがす

ケモインフォマティックス―予測と設計のための化学情報学
J.Gasteiger,T.Engel
丸善 / ￥ 18,900 ()
在庫あり。

SAR News No19に寄稿しました。この号はSVM,RFとかの統計手法を使ったSARから可視化とかグラフっぽい処理とか、ToolKitを使ったプログラミングとか、日本ではあまり見かけない内容なので面白いのではないでしょうか？

余談ですがハッカソンっぽい場もあるので、ちと書いてみるか（または飲んでみるか）という方がいれば連絡を下さい(まだ空きあります)。僕はOpenBabelを使ってGAMESSの量子化学計算で構造最適化をするラッパーを書く予定にしています。それかHaskellのchemoinformaticsライブラリのソース読んでます。

さて、今回の原稿を書くのにSphinxを使った。これ最高ですな、EPubも出力できるし、超便利。ただ、今回は原稿の提出がwordだったのでwordの行間調整したりとかよくわからない部分で難儀したけど(word嫌い)。でも履歴管理の仕組み（変更モードって言うの？）はメールでやり取りするときに便利だったけど。

あとバージョン管理はMercurialを使った。Gitに比べてwebアクセスが簡単に設定できたのでチョイス、あと文書管理だとそんなに複雑な操作しないし。

入門Mercurial Linux/Windows対応
藤原 克則
秀和システム / ￥ 2,310 ()
在庫あり。

以前、おもむろに思いついたLeadOptimizationのシミュレーターで、networkxはgml出力できてCytoscapeでimportすればいいじゃんということに気づいてヤル気が出た。

ついでにコードもちょっと直して、活性はpIC50とかそういうレンジにしてみて、μMオーダーから上げていくようにした。

エッジは世代を表していて、黄色から緑にむけてすすんでいく。ノードの色は活性を表していて黒から赤になるにしたがって活性が向上する。最初の方は活性が振れるけど段々向上していって安定するようにしてみた。

リードホッピングとは選択されなかったより良い未来を選択しなおすことであるとするならばそのような希望はLO初期のあたりに存在すると思うのだ。

以下、書いたコード

from random import random
import networkx as nx

class Branch(object):
    """LeadOptimization flow
    potency: pIC50 or such
    weight : range for activity
    """

    count = -1

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

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

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

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

if __name__ == "__main__":
    max_comps        = 500 # total compounds
    initial_branches = 1   # number of lead compounds
    lead_potency     = 5   # lead activity
    generation       = 0

    G=nx.Graph()

    heads = [Branch(lead_potency,2) for i in range(initial_branches)]
    map(lambda b: G.add_node(b.id, activity=b.activity),heads)

    while Branch.get_count() < max_comps:
        new_heads = []
        generation += 1
        for branch in heads[:int(2+5*random())]:
            for new_comp in branch.make_child(int(40*random()),branch.potency,branch.weight):
                G.add_node(new_comp.id, activity=new_comp.activity)
                G.add_edge(branch.id, new_comp.id, genneration=generation)
                if new_comp.activity > 7:
                    new_heads.append(new_comp)
        heads = new_heads
        heads.sort(key=lambda x:x.activity)
        heads.reverse()

    nx.write_gml(G,"test.gml")

MCSを求めるには比較したい二つの化合物のModular Productを求めて、それの最大クリーク探索をすればいい。

Modular Productとは

• 二つのグラフの属性の対応が同じノードが両方共隣接していてエッジの属性が同じ場合
• 二つのグラフの属性の対応が同じノードが両方共隣接していない場合

にエッジを張ったグラフ。

これを作ってあとはigraphの最大クリーク探索メソッドで

import openbabel as ob
from igraph import *

sm1 = 'OCC(N)=O'
sm2 = 'O=CC(N)=O'

obc = ob.OBConversion()
obc.SetInAndOutFormats('smi','smi')

source = ob.OBMol()      
obc.ReadString(source,sm1)
source_atoms = [a for a in ob.OBMolAtomIter(source)]

target = ob.OBMol()      
obc.ReadString(target,sm2)
target_atoms = [a for a in ob.OBMolAtomIter(target)]

pairs = []
for i in range(len(source_atoms)):
    for j in range(len(target_atoms)):
        if source_atoms[i].GetAtomicNum() == target_atoms[j].GetAtomicNum():
            pairs.append((source_atoms[i].GetIdx(),target_atoms[j].GetIdx()))

p = []
for x in range(len(pairs)-1):
    for y in range(x+1,len(pairs)):
        u = pairs[x]
        v = pairs[y]
        if u[0] != v[0] and u[1] != v[1]:
            sb = source.GetBond(u[0],v[0])
            tb = target.GetBond(u[1],v[1])
            if sb != None: sbo = sb.GetBondOrder()
            if tb != None: tbo = tb.GetBondOrder()

            if (sb == None and tb == None) or (sb != None and tb != None and sbo == tbo):
                p.append((pairs.index(u),pairs.index(v)))

g = Graph(p)
mc = g.largest_cliques()
for c in mc:
    print [pairs[i] for i in c]

実際にはline graph(エッジとノードを入れ替えたグラフ)のModular Productをもとめてクリーク探索したほうが計算量がかなり減るらしい。

が、今回はやっていない。

Handbook of Chemoinformatics Algorithms (Chapman & Hall/Crc Mathematical and Computational Biology)

Chapman & Hall / ￥ 8,895 ()
通常1～3週間以内に発送

Chemoinformatics and Computational Chemical Biology (Methods in Molecular Biology)

Humana Press / ￥ 14,386 ()
近日発売　予約可

MCSSTanimotoで組んだネットワークのメリットはなんといってもエッジの属性にMCSの情報(smilesとかInChi)をのっけられることであろう。

import sys
from Gaston import MCSSTanimoto,Gaston
import openbabel as ob

obc = ob.OBConversion()
obc.SetInAndOutFormats("sdf", "smi")

input = "pc_sample.sdf"
mol = ob.OBMol()
next = obc.ReadFile(mol,input)
mols = [mol]

while next:
    mol.StripSalts(14)
    mols.append(mol)
    mol = ob.OBMol()
    next = obc.Read(mol)

mols = mols[1:11]

for i in range(len(mols)-1):
    for j in range(i+1,len(mols)):
        sys.stdout.flush()
        mcs = MCSSTanimoto(mols[i],mols[j])
        if mcs.score() > 0.2:
            print "%s\t%s\t%s\t%s\t%2.3f" % (mols[i].GetTitle(), "mcs", \
            mols[j].GetTitle(), obc.WriteString(mcs.mcs).split()[0], mcs.score())

でも、実際に組んでみたら、解釈が意外に難しそうなことに気がついた。MCSSTanimotoのアルゴリズム自体が問題なのかもしれない。

MCSよりはmolblaster的に細切れにしてみて最頻出のフラグメントからネットワークを組んでいくか、ありがち置換基の定義済みテーブルを使って類似性の高い置換基動どうしをつなぐという経験ベースのアプローチのほうが直感的かもしれないなぁ。