狠狠撸

R+igraphではじめる
生物ネットワーク解析
竹本和広
九州工業大学情報工学研究院生命情報工学研究系
科学技術振興機構さきがけ
@kztakemoto
takemoto@bio.kyutech.ac.jp

第2回 Rでつなぐ次世代オミックス情報統合解析研究会@
Mar 8, 2013 理化学研究所横浜研究所

生物学の流れ
? 顕微鏡の発明による「細胞生物学」（17世紀）

? メンデルによる「遺伝学」（19世紀）

? X構造回折解析などによる「分子生物学」（20世紀）
? ワトソンとクリックによる「DNA＝遺伝子」

? Omics Science（20世紀後半）
? Genome, Transcriptome, Proteome, Metabolome, B
iome
? ヒトゲノム計画など

生体分子相互作用の
全体像が徐々に明らかに

異なる階層に、たくさんの役者

Biochemical Pathway Maps (http://web.expasy.org/pathways) (Zoology, 2nd edition, Brown Publishers, 1994)

データの蓄积

で、どうする？

Network Science
? Data-Driven Scienceの一種
? 統計解析?数理モデル
? 大量のデータを「知識」にする。

? 19世紀頃から社会学の一分野として認識されてい
た。
? Social network analysis
? Weak tieやStructural holeと富の関係など
? 最近ではSNS解析?サービスで利用

? 2000年頃から生物学で「精力的」に応用され始め
る。
? データが利用可能になったことにより

参考文献
? Barabási A-L, Oltvai ZN (2004) Network biology: understanding the
cell’s functional organization. Nat Rev Genet 5: 101–113.

? Barabási A-L, Gulbahce N, Loscalzo J (2011) Network medicine: a
network-based approach to human disease. Nat Rev Genet 12: 56–
68.

? Takemoto K (2012) Current understanding of the formation and
adaptation of metabolic systems based on network theory.
Metabolites 2: 429–457.

? Cho D-Y, Kim Y-A, Przytycka TM (2012) Chapter 5: Network biology
approach to complex diseases. PLoS Comput Biol 8: e1002820.

詳しくはWebの補足資料をご覧下さい
https://sites.google.com/site/kztakemoto/r-seminar-on-igraph---supplementary-information

このセミナーの目的

+igraphでネットワーク解析

? igraphの使い方の基礎

? 実際によく使われる手法
? 中心性解析
? ネットワークモチーフ
? コミュニティ抽出

http://igraph.sourceforge.net

なぜR+igraphか
? 利点
? Rユーザーは多い
? 他と比べてネットワーク解析の関数が豊富
? 他のデータと統合して解析できる。
? 欠点
? Visualizationは比較的劣る

? 他のネットワーク解析を支援するツール
? Pajek (Windowsのみで動作)
? NetworkX (Python)
? Cytoscape, Gephi (描画がきれい)

igraphパッケージの
インストールと読み込み

> install.packages("igraph")
いまからの説明はigraph ver. 0.6以上を想定しています。

バージョンを調べるには

> packageVersion("igraph")

パッケージの読み込み

> library(igraph)

ネットワークの読み込み

> d <- read.table("eco_EM+TCA.txt")
> g <- graph.data.frame(d,directed=T)
> g <- simplify(g,
+ remove.multiple=T,remove.loops=T

eco_EM+TCA.txt graph.data.frame:
Source Target データフレームからグラフに変換
2-Oxoglutarate Oxalosuccinate
2-Oxoglutarate Succinyl-CoA simplify:
2PG 3PG 多重エッジや自己ループの削除
2PG PEP
3PG 2PG
3PG BPG
AB
Acetyl-CoA Citrate …
Arbutin Arbutin_6-phosphate AB

他の形式からの読み込み

> g <- read.graph(filename,
+ format=“xxx”,
+ directed=T)

有名どころでは???
?edgelist (default) 01
?GML yEdやCytoscape 12
でもサポート
?GraphML 23
…
?Pajek

ネットワークを描画
> plot(g,vertex.label=V(g)$name,
+ vertex.size=5,
+ layout=layout.kamada.kawai)
V(g)$name:
ノード名の参照
layout:
レイアウトアルゴリズム
layout.fruchterman.reingold

layout.circle

layout.sphere
layout.grid

全体的なネットワーク特徴量を測る

? 次数分布

? クラスタ係数

? Assortative coefficient

? 最短パス

次数分布
? 次数（結合次数）：ノードが持つ枝数
次数kを持つノードの数
? 次数分布=
全体のノード数
> degree.distribution(g)
3
0.8 2
相対頻度 P(k)

0.6
0.4 2
0.2 1
0 2
1 2 3
結合次数 k

现実ネットワークの次数分布
スケールフリー性
アーキアバクテリア
? 次数分布がベキ関
数で近似できる。
? 明確な定義はない

? ランダムグラフと
は大きく異なる。

真核生物平均

Nature 407, 651 (2000) 例：代謝ネットワーク

次数分布とネットワークロバスト性

ハブが無いと弱い

ハブがあると偶発的
なアクシデントに強
い

しかし、攻撃に弱い

Scientific American 288, 50 (2003)

构造的な顽健性は生体システムの
頑健性と対応する場合がある
? 次数分布の均一性に注目
? 次数エントロピー

? 癌のシグナルネットワークの結合度の均一性と生
存率には負の相関 PNAS 109, 9209 (2012)

次数エントロピーの計算
> dist <- subset(degree.distribution(g),
+ degree.distribution(g)>0)
> -sum(dist*log(dist))

5年以内の生存率

クラスタ係数
? あるノードの近傍間においてエッジが張られ
る確率近傍ノード間
の辺数

平均クラスタ係数

> cc <- transitivity(g,type="localaverage",isolates="zero”)
> mean(cc)

Assortative coefficient
? 次数-次数相関の積率相関係数 Phys. Rev. Lett. 89, 208701 (2002)

r>0 r=0
> assortativity.degree(g)

r<0

r>0

r<0

フォールディング速度との関連
? コンタクトマップ
? ノード：アミノ酸

? エッジ：近接関係

Bioinformatics 23, 1760 (2007)

最短パス
a
距離行列 b
e
c
> average.path.length(g) d
aからdまでの最短経路を探す。
a b c d e
> sum <- a - 1 1 2 1
+ get.all.shortest.paths(g,
+ “a",”d",mode="out")
b 1 - 2 1 2
> V(g)[sum$res[[1]]] c 1 2 - 3 2
Vertex sequence:
[1]”a” “b” “d” d 2 1 3 - 1
> V(g)[sum$res[[2]]]
Vertex sequence:
e 1 2 2 1 -
[1]”a” “e” “d” > shortest.paths(g)

演习课题：代谢経路探索

D-Glucose から 2-Oxoglutarate
までのパスを見つけよう
> sum <- get.all.shortest.paths(g,
+ "D-Glucose","2-Oxoglutarate",mode="out")
> V(g)[sum$res[[1]]
Vertex sequence:
[1] "D-Glucose" "G6P” "F6P" "F16BP"
[5] "GAP" "BPG” "3PG” "2PG”...

> E(g)$color <- "grey"
> E(g,path=sum$res[[1]])$color <- "red"
> E(g,path=sum$res[[1]])$width <- 3

> plot(g,vertex.label=V(g)$name,
+ vertex.size=5,
+ layout=layout.fruchterman.reingold)

重要（そうな）なノードを見つける

? どのノードが重要で
しょうか？

? ネットワークにおける
位置と機能に関係は？

中心性解析
Soc Networks 1, 215 (1979)
Nature 411, 41 (2001)

良く知られた中心性（１）
? 次数中心性 > degree(g) / (vcount(g) - 1)
? 枝数の多い（次数の高い）頂点が重要

? 近接中心性 > closeness(g)
? その他の頂点に短い距離でつながる頂点が重要
頂点i-j間の最短距離

良く知られた中心性（２）
? 媒介中心性 > betweenness(g)
? 各ノード間の最短パスを考える時、よく通過する
頂点が重要
ノードj–k間の最短パスにおいてノードiを通る最短パスの数
ノードj–k間の最短パスの数

? 固有ベクトル中心性 > evcent(g)$vector
? 次数中心性の拡張版：近傍の中心性も考える

隣接行列つまり
頂点iの中心性

Google page rankの基本的な考え方

中心性の比较

Chpt. 2 in Statistical and Machine
Learning Approaches for Network
Analysis, John Wiley & Sons, 77 (2012).

Yeast PINにおいて相互作用数と
進化速度は負の相関がある
Science 296, 750 (2002)

相互作用の多いタンパク質の進化速度は小さい

PIN/GRNにおいて中心度の
高いタンパク質は生存に必須

PLoS Comput. Biol. 3, e59 (2007)

演習課題

(1) タンパク質相互作用ネットワーク
ecoli_ppi_Hu_etal_2009.txt
PLoS Biol 7, e96 (2009)

b3396 b2609
b4149 b3636
b4149 b3301

(2) タンパク質と
Nature 411, 41 (2001)
Essentialityの対応表
ecoli_proteins_essentiality_Baba2006MSB.txt
Mol Stst Biol 2, 2006.0008 (2006)

gene essential Essentialityと
Non-essential
b0001 N
… 中心性の関係に
Unknown
b0005 u
… ついて調べる
b0026 E Essential

次数中心性を用いる場合
> d <- read.table("ecoli_ppi_Hu_etal_2009.txt")
> g <- simplify(graph.data.frame(d,directed=F),
+ remove.multiple=T,remove.loops=T)

> ess <- read.table("ecoli_proteins_essentiality_Baba2006MSB.txt”,header=T)

> target <- intersect(ess$gene,V(g)$name)

> ess_select <- subset(ess,is.element(ess$gene,target)==T &
+ duplicated(ess$gene)==F)
> ess_select_ord <- ess_select[order(ess_select$gene),]

> deg <- degree(g)
> deg_select <- subset(deg,is.element(names(deg),target)==T)
> deg_select_ord <- deg_select[order(names(deg_select))]

b0001 E b0002 2 共通部分を
抜き出してソート b0001 E b0001 1
b0002 u b0004 5
b0002 u b0002 2
b0003 N b0001 1

> boxplot(log(deg_select_ord)~ess_select_ord$ess,ylab="log(Degree)")

结果
p < 2.2 × 10–16 using Wilcoxon test

Nature 411, 41 (2001)

Essential proteinは
比較的多くの相互作用を持つ

ネットワークモチーフ Science 298, 824 (2002)

実ネットワークにお
けるサブグラフの数

ヌルモデルにおけ
るサブグラフの数
の平均と標準偏差

>2なら有意に出現
(p=0.05)

ランダム化されたネットワーク
（ヌル仮説）に比べて頻出するパタン

ネットワークモチーフ
? フィードフォワードループなどが出現
? ネットワークのBuilding block
? 制御理論において、このモジュールは重要な
役割を果たす。

Science 298, 824 (2002)

ネットワークモチーフを
みつけるツール
? mfinder
Science 298, 824 (2002)

? mDraw
? by Uri Alon Group

MAVisto FANMOD 高速

Bioinformatics 21, 3572 (2005) Bioinformatics 22, 1152 (2006)

n-ノードサブグラフを数える
? 3ノードの場合（5以上は利用不可）
> c_real <- graph.motifs(g,3)
> c_real
[1] 0 0 2444 0 3014 10 180913 811 4 317
[11] 0 0 0 277 1 0

サブグラフの種類を確認する。 x in [0,15]（n=3の場合）
> plot(graph.isocreate(3,x),layout=layout.circle)

x=0 x=7 x=11 x=15

サブグラフを构成する
ノードの出力
? Feedforward loopの場合
> sum <- graph.get.subisomorphisms.vf2(g,graph.isocreate(3,7))
> sum
[[1]]
[1] 14 92 42

[[2]]
[1] 16 90 18
…

> V(g)[sum[[1]]+1] ノード名での出力
Vertex sequence:
[1] "caiF" "narL" "fnr"

Randomized Networks
? データ解析におけるヌルモデル
? あるネットワークの構造的特徴が有意なのかどう
かを検証する際に使われる。

? Configuration model
? Edge switching algorithm
? よく使われるが「適切な」ヌルモデルではないと
いう批判もある。
? 研究者の「適切さ」の数だけヌルモデル
? ヌルモデル戦争

Configuration Model Phys. Rev. E 71, 027103 (2005)

? 次数列を与える。
? Stubをランダムにつな
ぐ。
Chpt. 2 in Statistical and Machine
> deg <- degree(g) Learning Approaches for Network
> g_rand <- degree.sequence.game(deg) Analysis, John Wiley & Sons, 77 (2012).

連結性が保証されないことに注意
連結性を保証する場合にはvlモードを使う。
> deg_out <- degree(g,m="out")
> deg_in <- degree(g,m="in")
> g_rand <- degree.sequence.game(deg_out,deg_in,m=“vl”)

有向ネットワークの場合、相互エッジの数は保存されないことに注意

Edge Switching Algorithm
? 現実のネットワークを入力として任意に二辺を
選び、結合先を交換
> g_rand <- rewire(g,niter=n)
n: 辺を交換する回数
次数分布は変化しない
nはどのように設定すればよいのか？ Science 296, 910 (2002)
Science 298, 824 (2002)
終了条件「全てのエッジが少なくとも1回は交換されている」
≈クーポンコレクター問題
回、交換すると、全てエッジが少なくとも1回は交換されている確率

≈ e.g., c=5の時は0.993, c=10の時は0.99995

> rewire(g,niter=ecount(g)*(5+log(ecount(g)))/2)
有向ネットワークの場合、相互エッジの数は保存されないことに注意

演习课题：
ネットワークモチーフを見つける
大腸菌遺伝子制御ネットワーク
（RegulonDB ver 5.0）
regDB5.txt
実ネットワークのサブグラフを数える
Nucleic Acids Res 34, D394 (2006) graph.motifs(g,3)
acrR acrA
acrR acrB
acrR acrR Randomized networkを作る
ada ada
degree.sequence.game() or rewire()
Randomized networkのサブグラフを数える
注意

実際にはRandomized 平均と標準偏差からモチーフを見つけだす。
networkを作る際に、
相互エッジの数を保存
するという制約が必要

工夫の例
無向グラフに変換ランダム化して有向グラフに戻す

グラフを分ける

統合する
ランダム化

> g_unid <- delete.edges(g,subset(E(g),is.mutual(g,es=E(g))==T))
> g_bid <- delete.edges(g,subset(E(g),is.mutual(g,es=E(g))==F))
> g_bid_u <- as.undirected(g_bid)
> c_rand <- data.frame()
> for(i in 1:100){
+ g_rand_unid <- rewire(g_unid,niter=ecount(g_unid)*(50+log(ecount(g_unid)))/2)
+ g_rand_bid_u <- rewire(g_bid_u,niter=ecount(g_bid_u)*(50+log(ecount(g_bid_u)))/2)
+ g_rand_bid <- as.directed(g_rand_bid_u)
+ g_rand <- graph.union(g_rand_unid,g_rand_bid)
+ c_rand <- rbind(c_rand,graph.motifs(g_rand,3)) これだけでは不十分
+}
ほかにどんな工夫が必要でしょう

结果
> (c_real-colMeans(c_rand))/sqrt(colMeans(c_rand**2)-colMeans(c_rand)**2)
0 1 2 3 4 5 6
NaN NaN -6.8995100 NaN -2.8736176 -0.9491391 -3.0222612
7 8 9 10 11 12 13
6.8581263 0.4525696 -1.9818274 -1.3719887 -0.3144855 -0.4264014 2.4857428
14 15
1.6035675 NaN

> plot(graph.isocreate(3,7),layout=layout.circle)

> V(g)[sum[[1]]+1]
Vertex sequence:
[1] "caiF" "narL" "fnr”
[1] "chbR" "nagC" "crp“
[1] "csgD" "ompR" "crp”
[1] "csgD" "rcsA" "rcsB"

Feedforward loopが
ネットワークモチーフ

コミュニティ抽出
? 比較的密に連結した部分グラフ同士が疎に連
結しているような部分構造を見つけ出す。
? グラフクラスタリングに対応

どんなメリットがあるの？
? PPIや相関ネットワークにおい

て、機能的に類似な遺伝子が
同じコミュニティに属する
? 機能予測ができる。

? 複雑なネットワークの俯瞰的
な理解ができる

BMC Bioinformatics 7, 2 (2006)

大きく分けてふたつの形式がある

Overlapping

Non-overlapping
まずはこちらから紹介
Phys Rev E 74, 016110 (2006)

コミュニティ検出の基礎
? コミュ二ティ内の辺密度が高く、コミュティ
間の辺が疎であれば良い分割という前提
隣接行列（頂点 i と j が同じグルー
プに属していたら）
（そうでなければ）

グループのアサインを変え任意の次数列で構成されるラ
てQを最大化する。ンダムグラフの結合確率

Low Q High Q

コミュニティの最適な
アサインをみつける必要がある
? Qが最大となるc=(c1,c2,…cN)を見つける。
? 難しい（膨大な計算時間が必要）。

? 近似アルゴリズムを使う。

? 最適化アルゴリズムによってQが最大となるアサ
インを見つける。

アサイン変更
Q0<Q1: 採択
Q0>Q1: 棄却
Q0 Q1

有名な手法（１）
? 辺媒介性を用いた手法 PNAS 99, 7821 (2002)

? 媒介性の高い辺を削除して分割
> edge.betweenness.community(g) Phys. Rep. 486, 75 (2010)

? 貪欲アルゴリズムを用いた手法 Phys Rev E 70, 066111 (2004)

? Q値が高くなるように隣接（ノード?コミュニティ）
ペアをまとめていく。
> fastgreedy.community(g) (Q1–Q0) < (Q2–Q0)
Q1 Q2
Q0
採用

有名な手法（２）
Phys Rev E 74, 036104 (2006)
? スペクトル最適化法 PNAS 103, 8577 (2006)

? グラフラプラシアンを使ってQ値が最大となるよ
うな分割（カット）を探す。
> leading.eigenvector.community(g)

? 焼きなまし法（Simulated annealing）
Phys Rev E 74, 016110 (2006)
> spinglass.community(g)

アサインを
ちょっと変更 Q0<Q1: 採択
Q0>Q1:
の確率で採択
β（逆温度）を徐々に
Q0 Q1 大きくする（冷却）

計算结果の抽出
? 最大Q値
> max(data$modularity)
[1] 0.6141498

? 最大Q値におけるコミュニティ数
> max(data$membership)
[1] 5

? 最大Q値におけるコミュニティメンバーシッ
プ
> data$membership
[1] 4 1 1 5 2 2 1 5 2 1 2 3 2 2 1 1 4 3 5 4 5 5 2…

どれがいいの？
計算時間なら貪欲アルゴリズム
? 貪欲アルゴリズムが最小：n(ln n)2 Phys. Rep. 486, 75 (2010)

貪欲 < スペクトル法 < 辺媒介 < 焼きなまし法
（ネットワークが疎の場合）

どれがいいの？
パフォーマンスなら焼きなまし法

辺媒介 < 貪欲 < スペクトル法 < 焼きなまし法
Phys. Rep. 486, 75 (2010)

Q値の問題点
? Resolution limit 解像度限界 PNAS 104, 36 (2007)

? 小さなコミュニティを見逃してしまうことがあ
る。
? 疎なネットワークな場合がおおよそこれに当てはま
る。 m=5
m-ノード完全グラフ m=3
? コミュニティ内の次数の合計Kが

m>4のとき、それぞれの完全グラフ
がコミュニティとして認識される。

モジュラリティのより良い定義が必要

より良い指標：
Modularity Density Phys Rev E 77, 036109 (2008)

? Simple definition: コミュニティ外の平均次数

コミュニティ内の平均次数
D1
? More general definition

<
where

D2

演习课题：
代謝経路のコミュニティ抽出
eco_EM+TCA.txt

edge.betweenness.community(g)
fastgreedy.community(g)
leading.eigenvector.community(g)
spinglass.community(g)

貪欲アルゴリズムを使う場合
?
d <- read.table("eco_EM+TCA.txt")
g <- simplify(
+ graph.data.frame(d,directed=F),
+ remove.multiple=T,remove.loops=T)
data <- fastgreedy.community(g)
V(g)$color <- data$membership

结果
> plot(g,vertex.size=5,
+ vertex.label=V(g)$name, デンドログラムの表示
+ layout=layout.kamada.kawai) > dend <- as.dendrogram(data)
> plot(dend)
解糖系?糖新生の上部

=>TCA
解糖系?糖新生下部

TCA (酸化)
TCA (還元)

Overlapping Communities
? コミュニティがオーバーラップしている（階
層性を持っている）と考える方がより自然

? Clique percolation
Nature 435, 814 (2005)

? Edge partition
Phys Rev E 80, 016105 (2009)

? Link communities
Nature 466, 761 (2010)

? Overlapping cluster generator
Nature 435, 814 (2005) Bioinformatics 28, 84 (2012)
利用可能なのはふたつ

> library(linkcomm)

Clique Percolation Method Nature 435, 814 (2005)

? k-clique隣接に基づいたパーコレーションか
らコミュニティを同定する。
? CFinderで計算できます (cfinder.org)
k-clique隣接
Yes No

k=3の場合
k=4の場合

Edge Partition Phys Rev E 80, 016105 (2009)

? Line graphに変換してからコミュニティ抽出
Line graph変換

Softwareは利用可能です
https://sites.google.com/site/linegraphs/

エッジのクラスタリング

Link Communities Nature 466, 761 (2010)

（１）エッジ間の類似性計算（２）適当な階層的クラスタリング法
でクラスタリングする

k
eik ejk
i j

（３）分割密度Dが最大になる分割を探す
（コミュニティ内が密になるように）

論文では単連結法が使用されている

> data <- getLinkCommunities(g)

Link Community Method Shows High
Composite Performance
Manually-curated Clustering:
研究者（コンソーシアム）が決定した分類
例：Gene OntologyやKEGG Pathway

Overlapping Cluster Generator

1. 適当な初期状態を設定

2. Q値（オーバーラップ
版）が高くなるように
オーバーラップしている
コミュニティを融合する
（貪欲アルゴリズム）

CFinderやLinkCommよりも
良いパフォーマンスを示
す。
> data <- getOCG.clusters(g)

演习课题：代謝経路のコミュニ
ティ抽出（Overlapを考慮）
eco_EM+TCA.txt

linkcommはネットワークの読み込みが異なる
?
> g <- read.table("eco_EM+TCA.txt")

> data <- getOCG.clusters(g)

> plot(data,type="graph”)
> data$clusters

结果：LinkCommの場合
任意の分割密度での结果を表示
> lc <- newLinkCommsAt(data,cutat=0.84)
> plot(lc,type="graph")

1本のエッジで構
成されるコミュ
> plot(data,type="graph")
ニティは表示さ
れない

结果：OCGの場合
糖の取り込み

> plot(data,type="members")

解糖系?
糖新生
上部 TCA（還元） TCA
（酸化）

複数のコミュニティに
属するものだけ表示

解糖系?糖新生下部
=>TCA

まとめ
? 中心性解析
? ネットワークと生体分子機能の関係を調べる。

? ネットワークモチーフ
? 意図的に構成された部分構造をみつける。

? コミュニティ抽出
? 機能予測や俯瞰的な理解を助ける。

@kztakemoto
takemoto@bio.kyutech.ac.jp

狠狠撸

R seminar on igraph

More Related Content

R seminar on igraph