Drkcore

構造最適化シミュレーションをCytoscapeで

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

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

projectsim

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

リードホッピングとは選択されなかったより良い未来を選択しなおすことであるとするならばそのような希望は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")

2010-09-03T20:10:29

igraphを使ってmaximum common substructure(MCS)を求める

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 ()
近日発売　予約可

2010-08-18T18:44:45

LeadOptimizationをシミュレーションしたい

ふと、つぶやいた。

例えば

プロジェクトで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")

project

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

2010-03-10T21:56:46

禁断探索法をpythonで

ヒューリスティクスは割と好きだ。スパッととけないものをなんとか現実と折り合いを付けながらよりよい解を見つけるみたいな、なんとも煮え切らない美しさのなさがたまらん。インフォマティクスも似たような香りはする。

特に、「シミュレーションできない学問は学問として未成熟である」という言葉に従うのであれば、創薬研究は（学際領域）と言われている割には未成熟な学問の組み合わせのために発見的な手法の割合が増えすぎるし、精度の高い予測法の誕生もまだまだ先であろう。というわけで、発見的な探索手法は当分有用だし、学際領域故に融合部分での応用が期待されるため、アルゴリズムとしてきちんと押さえておくと色々役に立つ。

メタヒューリスティクスの数理の4章はpythonのコードが載っているので、勉強になる。

メタヒューリスティクスの数理
久保幹雄,J. P. ペドロソ
共立出版 / ￥ 3,675 ()
在庫あり。

第4章応用
4.1 グラフ分割問題
4.2 最大安定集合問題
4.3 グラフ彩色問題
4.4 巡回セールスマン問題
4.5 2次割当問題
4.6 多制約ナップサック問題
4.7 数分割問題

4.1のグラフ分割問題をnetworkxで可視化するように書き換えた。

path

import math,random
from sys import maxint as Infinity

def rnd_graph_fast(nnodes, prob):
    nodes = range(nnodes)
    adj = [[] for i in nodes]
    i=1;j=1
    logp = math.log(1.0 - prob)
    while i<nnodes:
        logr = math.log(1.0 - random.random())
        j += 1+int(logr/logp)
        while j>=i and i<nnodes:
            j -= i
            i += 1
        if i<nnodes:
            adj[i].append(j)
            adj[j].append(i)
    return nodes, adj

def construct(nodes):
    sol = [0 for i in nodes]
    for i in range(len(nodes)/2):
        sol[i] = 1
    return sol

def evaluate(nodes, adj, sol):
    cost = 0
    s = [0 for i in nodes]
    d = [0 for i in nodes]
    for i in nodes:
        for j in adj[i]:
            if sol[i] == sol[j]:
                s[i] += 1
            else:
                d[i] += 1
    for i in nodes:
        cost += d[i]
    return cost/2, s, d\

def find_move(part, nodes, adj, sol, s, d, tabu, tabulen, iteration):
    mindelta = Infinity
    istar  = None
    for i in nodes:
        if sol[i] != part:
            if tabu[i] <= iteration:
                delta = s[i] - d[i]
                if delta < mindelta:
                    mindelta = delta 
                    istar = i
    if istar != None:
        return istar, mindelta
    tabu = [0 for i in nodes]
    return find_move(part, nodes , adj, sol, s, d, tabu, tabulen, iteration) 

def move(part, nodes, adj, sol, s, d, tabu, tabulen, iteration):
    i, delta = find_move(part, nodes, adj, sol, s, d, tabu, tabulen, iteration)
    sol[i] = part
    tabu[i] = iteration + tabulen
    s[i], d[i] = d[i], s[i]
    for j in adj[i]:
        if sol[j] != part:
            s[j] -= 1
            d[j] += 1
        else:
            s[j] += 1
            d[j] -= 1
    return delta

def tabu_search(nodes, adj, sol, max_iter, tabulen):
    cost, s, d = evaluate(nodes, adj, sol)
    tabu = [0 for i in nodes]
    bestcost = Infinity
    for it in range(max_iter):
        cost += move(1, nodes, adj, sol, s, d, tabu, tabulen, it)
        cost += move(0, nodes, adj, sol, s, d, tabu, tabulen, it)
        if cost < bestcost:
            bestcost = cost
            bestsol = list(sol)
    return bestsol, bestcost


if __name__ == "__main__":
    import networkx as nx
    import matplotlib.pyplot as plt
    num_nodes = 50
    nodes, adj = rnd_graph_fast(num_nodes,0.3)
    G=nx.Graph()
    for i in range(num_nodes):
        G.add_node(i)
    for i in range(num_nodes-1):
        for j in adj[i]:
            if i < j:
                G.add_edge(i,j)
    pos=nx.spring_layout(G)

    max_iter = 1000
    tabulen = 10
    sol = construct(nodes)
    sol, cost = tabu_search(nodes, adj, sol, max_iter, tabulen)

    rnodelist = []
    bnodelist = []
    for i in range(len(sol)):
        if sol[i] == 1:
            rnodelist.append(i)
        else:
            bnodelist.append(i)
    nx.draw(G,pos,nodelist=rnodelist,node_color='r')
    nx.draw(G,pos,nodelist=bnodelist,node_color='b')
    plt.savefig("path.png")

これをちょっと変えてSimulated Annealingにできる。

2010-01-12T10:34:45

macbookにgraphviz,pygraphviz,networkxを入れた

コンパートメントモデルを書くのにいちいちInkspaceを立ち上げるのもめんどうなのでGraphvizにでも頼るかと思ったらmacbookには入ってなかった。networkxもついでに入れといた。

sudo port install graphviz
sudo easy_install-2.6 pygraphviz
sudo easy_install-2.6 networkx

使う

>>> import pygraphviz as pgv
>>> d={'1': {'2': None}, '2': {'1': None, '3': None}, '3': {'2': None}}
>>> A=pgv.AGraph(d)
>>> A.layout(prog='dot')
>>> A.draw('pyg.png')

test