cytoscape / Drkcore

Mishima.syk #4やりました

Sun, 26 Oct 2014 16:38:07 +0919

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

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

反省点

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

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

私のスライド

Mishimasyk141025 from Kazufumi Ohkawa

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

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

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

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

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

週末の僕のハンズオンはCytoscape, iGraph, RDKitあたり

Tue, 21 Oct 2014 20:47:23 +0919

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

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

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")

今週末のMishima.sykは化合物ネットワークが多めになりそう

Mon, 20 Oct 2014 21:56:17 +0919

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

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

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

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

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

CytoscapeでDate型を扱う

Tue, 12 Aug 2014 23:35:42 +0919

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

要するにunixtimeですね。

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

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

Cytoscapeでchemoinformatics

Sun, 10 Aug 2014 12:29:58 +0919

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)])))

散布図を密度の情報からグラフに変換する

Sat, 04 Aug 2012 13:03:53 +0919

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

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

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

アルゴリズム

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

例

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

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

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

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

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

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

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

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

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

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

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

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

選曲ネットワークをCytoscapeで

Mon, 23 May 2011 20:22:58 +0919

ゴールデンウィークを利用して250ノード、220エッジくらいまで選曲データを溜めたので、Cytoscapeで描画してみた。SIFフォーマットにして読み込んでみた。

全体像はだいたいこんな感じ。たまにacid jazzとかかけてたので、それが孤立したネットワークとして左下に表示されてるが、基本はドラムンベース。

中心部を拡大してみた。

sebaはかなり好きななんだが、Audioもよく聴いている。

ところで、文字化けしてんのはどうやったらいいんだろうか？

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

Sat, 09 Oct 2010 10:14:57 +0919

以前、おもむろに思いついた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で組んだネットワークのメリット

Mon, 16 Aug 2010 20:55:30 +0919

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

引き続き化合物のネットワークをいじっている

Sat, 14 Aug 2010 09:13:34 +0919

ノード間で類似性が最大となるような関係にひとつだけエッジを張るようにしてみた。

#!/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.