ゴールデンウィークを利用して250ノード、220エッジくらいまで選曲データを溜めたので、Cytoscapeで描画してみた。SIFフォーマットにして読み込んでみた。
全体像はだいたいこんな感じ。たまにacid jazzとかかけてたので、それが孤立したネットワークとして左下に表示されてるが、基本はドラムンベース。
中心部を拡大してみた。
sebaはかなり好きななんだが、Audioもよく聴いている。
ところで、文字化けしてんのはどうやったらいいんだろうか？
以前、おもむろに思いついた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")
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的に細切れにしてみて最頻出のフラグメントからネットワークを組んでいくか、ありがち置換基の定義済みテーブルを使って類似性の高い置換基動どうしをつなぐという経験ベースのアプローチのほうが直感的かもしれないなぁ。
ノード間で類似性が最大となるような関係にひとつだけエッジを張るようにしてみた。
#!/usr/bin/env python
# -*- encoding:utf-8 -*-
import pybel
mols = list(pybel.readfile("sdf", "pc_sample.sdf"))
fps = [x.calcfp() for x in mols]
for i in range(len(fps)-1):
max_sim = 0
max_name = ""
for j in range(i+1,len(fps)):
sim = fps[i] | fps[j]
if sim > max_sim:
max_sim = sim
max_name = mols[j].title
print "%s\t%s\t%s\t%2.3f" % (mols[i].title, 'sim', max_name, max_sim)
ついでに類似度属性で色が変わるようにしてみた。
わからん。さっきのエッジに加えて類似性が高い化合物同士のエッジは加えてみた。
for i in range(len(fps)-1):
max_sim = 0
max_name = ""
targets = []
for j in range(i+1,len(fps)):
sim = fps[i] | fps[j]
if sim > max_sim:
max_sim = sim
max_name = mols[j].title
if sim > 0.7:
print "%s\t%s\t%s\t%2.3f" % (mols[i].title, 'sim', mols[j].title, sim)
if len(targets) == 0:
print "%s\t%s\t%s\t%2.3f" % (mols[i].title, 'sim', max_name, max_sim)
ちょっとまとまりが出てきたけど、あれだなぁ。もとのsdfがそもそもあかんのかなぁ。
で、最小全域木(MST)で描くことも考えたけど化合物間の類似性の最小経路を求めてもなぁ、、、とか思ってヤル気がおきない。類似性でつないだネットワークのMSTの意義ってなんなんだろうなぁと。
なんやろかー、なんやろかー。論文もう一回ちゃんと読んでみようかなぁ。
chemoinformaticsでもネットワークアナリシスは重要だと思う。今はマイナーだけど今後精力的に研究されんじゃないかなぁと。
pybel使うとシミラリティベースのネットワークは簡単に作れる。
import pybel
mols = list(pybel.readfile("sdf", "pc_sample.sdf"))
fps = [x.calcfp() for x in mols]
for i in range(len(fps)):
ids = []
for j in range(i,len(fps)):
if (fps[i] | fps[j]) > 0.5 and i != j:
ids.append(mols[j].title)
print "%s %s %s" % (mols[i].title, "sim", " ".join(ids))
SIF(simple interaction format)形式なので、拡張子sifでcytoscapeに読み込まさせればOK。
あとでケモインフォクックブックにも載せておこう。
pubchemはクラスタリングの結果をgml形式でダウンロードできるので、それをcytoscapeで読み込んで解析することができる。早速データをとってくる。benzothioleで検索すると大体3000件くらいヒットするのでこれをクラスタリングしたものをgmlでダウンロードして、cytoscapeで読み込んでみた。
適当に描かせて眺めてみる。なんかネットワークっぽくなってますの。
さらに、ズームしてみる。
うーん、クラスタリングの結果を見てもちょっと面白くない。むしろ部分構造(MCS：Maximum Common Substructure)から構築したネットワークのほうがわかりやすいだろうな。
あとは、合成の時系列でネットワークをつくればいいと思うが、その場合、元々参考にした化合物の構造への情報をケミストから聞かなきゃいけないのでちょっとめんどいなぁ。
ということはpatentか？あれは先行技術の特許番号ついてたはずだからあれでネットワーク構築すればいいのか。これだったら、描いてみておもしろそう。
TODO: 奇麗な絵が描けるようにcydocなどを参考にちょっと勉強しとく