参加されたみなさんお疲れ様でした。

僕はNGSの解析の話が聞けたし、ネットワークまわりの色々なヒントが聞けて非常に勉強になりました

反省点

• 懇親会は予約する (脱香香飯店w)
• 告知は余裕を持ってする
• 反省会はちゃんと反省をする

反省会は忘年会を兼ねてやればいいですかね（やごみで）…ｗ

私のスライド

自分のスライドというかハンズオンのセッションは、Cytoscape以外に色々いれないといけない感じになってしまったので、予めサンプルを用意しておけばよかったかなぁと思った。

これからigraphを使ったネットワーク分析をしてみようと思う人が増えてもらえると嬉しいなと。

ネットワーク分析 (Rで学ぶデータサイエンス 8)
鈴木 努
共立出版 / 3564円 ( 2009-09-25 )

懇親会は(安定の？)香香飯店

二次会はワインバーみたいなとこ（名前忘れた）

RDKItはインストールが面倒くさいのと、実務だとPipelinePIlotとかOpenEye使うだろうからさらっと流すかも。その分iGraphのいじり方を多めに説明するかもしれない。

まだスライドに手をつけてないからわからないけどｗ

• Mishima.syk #4
• Mishima.syk #4 懇親会

iGraphを使えばノードとエッジに属性を簡単に付加できてgml形式で出力できる。このファイルはすぐにCytoscapeで開くことが出来て便利です。sifとかタブ区切りテキストは色々と面倒ですね。

from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs, MolToSmiles
import numpy as np
from igraph import Graph
import numpy as np
import logging
logging.basicConfig(level=logging.DEBUG)

threshold = 0.5

suppl = Chem.SDMolSupplier('syk.sdf')

fps = [AllChem.GetMorganFingerprintAsBitVect(mol, 2) for mol in suppl]
smiles = [MolToSmiles(mol) for mol in suppl]

mat = []
for i in range(len(fps)):
    mat.append([DataStructs.FingerprintSimilarity(fps[i],fps[j]) for j in range(len(fps))])

sm = np.array(mat)

def search_root(sm, remain):
    cand = None
    max_c = 0
    for l in remain:
        num_c = len(np.where(sm[l] > threshold)[0].tolist())
        if num_c > max_c:
            max_c = num_c
            cand = l
    return cand

def check_edge(sm, current, root):
    cand = None
    max_c = 0
    candidates = set(np.where(sm[root] > threshold)[0].tolist())
    for c in candidates:
        if c in current:
            if sm[root, c] > max_c:
                max_c = sm[root, c]
                cand = c
    return cand

edges = []
sim_edges = []
np.fill_diagonal(sm, 0)
remain = set(range(len(sm[0])))
current = set()
new = set()
ith = 1

while remain:
    logging.debug("### {} th ###".format(ith))
    ith += 1
    root = search_root(sm, remain)
    if root is None:
        break
    logging.debug("root: {}".format(root))

    cand = check_edge(sm, current, root)
    if cand:
        logging.debug("connect {} and {}".format(cand, root))
        edges.append((cand, root))
        current.remove(cand)
        sm[:, cand] = 0

    remain.remove(root)
    sm[:, root] = 0

    for i in np.where(sm[root] > threshold)[0].tolist():
        if i in remain:
            edges.append((root, i))
            sim_edges.append(sm[root, i])
            new.add(i)

    logging.debug("new nodes: {}".format(new))
    current |= new
    remain -= current
    new = set()

g = Graph(edges)
g.vs["smiles"] = smiles
g.es["similarity"] = sim_edges
g.save("test.gml")

今回はなかなか企業にいると触れることのないと噂のNGSデータの分析のハンズオン(@bonohu)が目玉ですが、他には化合物のネットワーク関連の話が多そうなので興味があれば参加してください。懇親会も、ディープ目のところを攻めていくので出ると楽しいハズｗ

• Mishima.syk #4
• Mishima.syk #4 懇親会

尚、私の発表は、Cytoscape.jsで化合物ネットワーク解析の触りをやろうというハンズオンにする予定です。ハンズオンが終わったらこんなのを描けるようになるのが目標です。

とかいいつつ、エッジの張り方を色々と試行錯誤しはじめていた。下はある構造最適化の論文の化合物群に対して最近傍同士をつないでいったネットワーク。

なんか細切れになっていて気持よくない。やっぱルートから一本につながってないとあかんよなとｗ

どういうネットワークを作るかというのが非常に重要だと思っているので、距離に使う記述子の話とか、どういう風にノード同士をつなげるかとかそういう話もしてみたいかな。

時系列データを扱うときには日付で色の濃度を変えたりしたいことが多いけどCytoscapeにはDate型が用意されてない？のでintに変換して代用しようとした。

要するにunixtimeですね。

int(time.mktime(time.strptime("2014-08-12 18:00:00", "%Y-%m-%d %H:%M:%S")))

とやればint型に変換されるのでめでたしめでたしと言いたいところだがいつの日付かさっぱりわからないｗ

Chemviz2を使い始めたのでメモ

インストール

• Cytoscapeのサイトから3.xの最新バージョンをダウンロードしてインストール
• Apps -> App ManagerからChemviz2を探してインストール

使ってみる

pychembldbを使ってノード用の属性ファイルとnetworkファイルを作った。ネットワークはとりあえずランダムにつないでみた。

Cytoscapeを起動して"network.sif"をインポート。続いてテーブルファイルとして"node.csv"をimport プライマリキーはsmilesにする。

構造描画の設定はApps -> Chemoinformatics Tools -> SettingsでSmiles Attributesをnode.shared.nameを選んでおく。 構造描画するときには右クリックしてオプションダイアログのApps -> Cheminformatics Toolsでpaint allかselectを選ぶ。

サイズをpIC50にして、オーガニックレイアウトで表示してみた。なかなかいい感じ

よく分からなかったところ

• nameとshared nameの違いがよく分からなくてドキュメントあたってみたけど結局良くわからなかった
• Chemviz2でsmilesをshared nameにしなきゃいけなくてsifファイルをsmilesで作ったのがダサかった
• (sifのidはchembl_idとかにしたかったけど)

コード

from pychembldb import *
from copy import copy
from random import shuffle, random
from math import log10
#Inhibition of recombinant Syk
#Bioorg. Med. Chem. Lett. (2009) 19:1944-1949
assay = chembldb.query(Assay).filter_by(chembl_id="CHEMBL1022010").one()

with open("node.csv","w") as f:
    f.write('"ID","pIC50","ALOGP","MWT","SMILES"\n')

    for act in assay.activities:
        if act.standard_relation == "=":
            f.write('"{}",{},{},{},{}\n'.format(act.compound.molecule.chembl_id,
                                            9 - log10(act.standard_value),
                                            act.compound.molecule.property.alogp, 
                                            act.compound.molecule.property.mw_freebase,
                                            act.compound.molecule.structure.canonical_smiles))

with open("network.sif", "w") as f:
    smiles_list = [ act.compound.molecule.structure.canonical_smiles for act in assay.activities if act.standard_relation == "="]
    for i in range(len(smiles_list)):
        tl = copy(smiles_list)
        tl[i:i+1] = []
        shuffle(tl)
        f.write('{} randomnetwork {}\n'.format(smiles_list[i], " ".join(tl[:int(random() * 5)])))

複数のパラメータが相互作用しあっている状態で最適なパラメーターセットを探索する際に、あまり精度の高いフィードバックが得られないと最適解の近傍の密度が高くなり、さらに準最適解があちこちに散らばるようなことがおこる。

それはPCAとかICAで空間を適当にとってあげることが出来て、特に第一主成分と第二主成分を平面にとってプロットしてやると、散布図として表現できる。

さてこれを、グラフにしてやればCytoscapeで扱えて便利だろうと思ったことがあったのでやり方を考えてみた。

アルゴリズム

1. 任意の点から非類似度Tの距離にある点の数をすべての点について数える
2. 最も点を含むものが多い点をrootとし、含まれる点との間で親子関係をつくる
3. 2の点を除いて1を行い新たな親子を決める。その際に親から非類似度T内にある既存の子のうち一番近いものを結ぶ（階層ができる）
4. これをすべての点がなくなるまで繰り返す。

例

以下の様な散布図を考える。

総当りで近傍の数をかぞえると、次の5つの点が含まれる空間が一番密度が高いことが分かる。

中心と含まれる点に対して親子関係をつくる(P1, C1-C4)。

これらを除いて、近傍の数をかぞえると次に密度が高いのが右側で3つの点を含んでいる。

これらに関しても先ほどと同様に親子関係をつくっておくが、新しくできた親に関しては既存の子とつなげられないかをチェックする。

するとC2 -> P2をつなげるのが良いことがわかる。

同様に行う（図の左側）。

これもC4 -> P3とつなげることができる。

最終的に全部つなげることができた。

実際に階層型の木を描くと以下のようになる。

単純だけど、全体の構造を保ったままCytoscapeで眺められていいかなーと思うんだけど。

もっといい方法あったら教えてもらえると嬉しいです。

ゴールデンウィークを利用して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の意義ってなんなんだろうなぁと。

なんやろかー、なんやろかー。論文もう一回ちゃんと読んでみようかなぁ。

• Finding key members in compound libraries by analyzing networks of molecules assembled by structural similarity.