R+igraph

Tutorial on R+igraph - supplementary information

This document was originally prepared for a Japanese workshop on the software R; however, it has been updated through interaction.

Especially, part of the tutorial is updated by Keiichiro Ono (Cytoscape Core Developer at UCSD) for cyREST, an app for driving Cytoscape from R/Python/Julia/Node.js/etc. The tutorial, updated by Ono, is very useful for people who mainly use English and/or are interested in Cytoscape+R+igraph.

if you have any questions or requests, please feel free to contact me via Twitter and e-mail.

はじめに

内容や使い方に疑問・質問等ありましたら。Twitterやメールでお気軽にご連絡下さい。

第2回Rでつなぐ次世代オミックス情報統合解析研究会「R+igraphではじめる生物ネットワーク解析」の補足情報です。 したが、その後の交流で色々と改訂され、igraphの使い方のページになっています。記載される内容は以下のようです。

• ネットワークを読み込む

• ネットワークを描画する

• ネットワークにおける各ノードの指標を計算する：中心性（Centrality）解析など

• ネットワークの全体的な特徴を計算する：クラスタ係数、平均パス長、経路検索など

• ネットワークのランダム化（ヌルモデル）：帰無仮説の作成

• コミュニティ検出（Community detection）

  • コミュニティが重複しない場合（Non-overlapping community detection）

  • コミュニティが重複する場合（Overlapping community detection）

  • Functional cartography：コミュニティ情報を用いてノードの役割を分類する

• ネットワークモチーフ：ネットワークに有意に出現するサブグラフを見つける

• ネットワークを比較する

• その他応用

  • Box-covering algorithm：ネットワークのフラクタル次元を測る

  • コア－周辺構造、蝶ネクタイ（入力－コア－出力）構造の検出

  • Threshold selection：ランダム行列理論に基づく閾値選択で相関ネットワークを構築

  • 最小支配集合を求める：ネットワーク制御可能性との関連

  • 最大マッチングを求める：ネットワーク制御可能性との関連

Rパッケージのひとつであるigraphではネットワーク解析に必要なたくさんの関数が登録されており、ネットワーク解析を簡単に行うことができます。ネットワーク解析を支援するツールはたくさんありますが、igraphは多くの解析手法が利用可能であることに加え、Rとの併用で、ネットワーク解析の結果を他のデータと統合した柔軟な解析ができるという利点があります。

• 演習で使用するスクリプトとデータはこちらからダウンロードして下さい。

• 発表スライドはこちらからご覧頂けます。

【演習問題】CodeIQに「友好関係ネットワークから派閥を検出」と題して、コミュニティ抽出に関する問題を提供しました。演習問題として、自身の理解の確認として取り組んではいかがでしょうか。解答の受付は既に締め切られていますが、ソースコードをお送り頂けばフィードバック致します。なお、CodeIQ Blogでこの問題の解説を行っております（解答が記載されていますので、ご注意下さい）。ここではネットワーク分析に関する周辺知識についても記載しておりますので、より深い理解に役立つと思います。

以下の内容は@sato_mitsunori様のブログでも紹介されています。重要な点が簡潔にまとめられ、図も豊富でとても分かりやすい良い資料です。こちらも参考にして頂ければ幸いです。

以下の内容の一部はCytoscapeの開発者である大野圭一朗さん（UCSD）によって、cyREST（CytoscapeとR/Python/Julia/Node.jsなどとの連携のためのアプリ）のチュートリアルとして、改訂されています。Cytoscapeとの連携に興味がある方は【こちら】をご参考ください。

この演習ではigraphのバージョンが0.6以上であることを想定しています（igraphのバージョンが0.6未満の場合、スクリプトの動作は保証されません）。以下のようにしてigraphのバージョンを確認してください。

> packageVersion("igraph")

バージョンが0.6未満、もしくはインストールされていない場合は以下のようにして最新版をインストールして下さい。

• Rコンソールを立ち上げる。

• メニューバーから「パッケージとデータ」＞「パッケージインストーラ」を立ち上げる

• 「パッケージ検索」にigraphと入力して、「一覧を取得」を押す。

• igraphを選択して、「選択をインストール」を押す。

もしくはターミナル上でRを起動させ、以下のコマンドを実行して下さい。

> install.packages("igraph")

部分的にlinkcommというパッケージも使用しますので、これも併せてインストールして下さい。

> install.packages("linkcomm")

linkcommは以下のパッケージも要求します。これらについてもインストールしてください。

> install.packages("igraph0") # 必要なくなりました。

> install.packages("RColorBrewer")

生物ネットワーク解析（大きく言って、Network Biology）に関しては以下の文献が入門になります（公開されている／オープンアクセスである論文のみリストアップしています）。

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

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

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

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

igraphの基本的な使い方

igraphの基本的な使い方をまとめています。演習で全てについて言及する訳ではありませんが、一度目を通しておくと、理解がより深まると思います。

各関数のオプションなどの詳細についてはigraphのドキュメントも参照してください。

igraphライブラリの読み込み

> library(igraph)

ネットワーク（グラフ）データの読み込み

以下のように、ノード番号が0...Nと割り振られているエッジリストがあるとします。

0 1

1 2

...

3 N

このような場合、以下の様にしてネットワークを読み込むことができます。オプションdirectedはエッジの方向性を意味します。 directed=Fで、方向性の無いエッジで構成されるネットワーク（無向グラフ）として読み込まれます。

> g <- read.graph(filename,directed=F)

また、read.graph関数はGML: Graph Modeling LanguageやGraphML（グラフ専用のXMLベースのファイルフォーマット）といった有名ないくつかのフォーマットでの入力が可能です。GMLの場合はオプションformatを用いて以下のように行います。ファイルに、エッジの方向性の有無（無向／有向きグラフ）やエッジの重み、ノードのラベルなどの指定がある場合は、それらも自動的に読み込まれます。

> g <- read.graph(filename, format="gml")

GraphMLの場合はgraphmlを指定することで、読み込みが可能です。他にも、Pajek（pajek）などが使えます。

より現実的には、以下のようにノード名が文字列で与えられているエッジリストを入力とする場合が多いでしょう。1列目と2列目が二項関係を示し、3列目が重みだとします。

A B 1

B C 3

D C 4

...

このような場合は以下のようにして読み込むことができます。

> d <- read.table(filename)

> g <- graph.data.frame(d[1:2],directed=F)

エッジの重みは以下のようにして読み込む。

> E(g)$weight <- d[[3]]

以下のように行列データが与えられている場合を考えます。

a b c d

a 1 3 0 0

b 1 0 2 0

c 0 1 0 1

d 0 0 4 0

このような場合は以下のようにして読み込むことができます。

> d <- as.matrix(read.table(filename,header=T,row.names=1))

# 重み付き有向ネットワークの場合

> g <- graph.adjacency(d,mode="directed",weighted=T)

# 無向ネットワークの場合

> g <- graph.adjacency(d,mode="undirected",weighted=NULL)

ノードのラベルは以下のようにして参照できる

> V(g)$name

ノードラベルによるエッジリストの場合、エッジを持たないノードが無視されてしまいますが、ノードラベルの対応表を用意しておけば大丈夫。

> ver <- read.table(filename)

> g <- graph.data.frame(d,directed=F,vertices=ver)

ノード間の重複したエッジや自己ループ（自分から自分につながるエッジ）を除く

> g <- simplify(g0,remove.multiple=T,remove.loops=T)

まとめて書くと便利

> g <- simplify(read.graph(commandArgs()[6],directed=F),remove.multiple=T,remove.loops=T)

> g <- simplify(graph.data.frame(d,directed=F),remove.multiple=T,remove.loops=T)

ネットワークの描画

> plot(g)

オプションを付ければより高度な描画が可能になります。以下の様に様々なオプションが利用可能です。

> plot(g,

+ vertex.size=15, #ノードの大きさ

+ vertex.shape="rectangle", #ノードの形

+ vertex.label=V(g)$name, #ノード属性nameをノードラベルにする。

#ノード属性Factionを用いてノードに色づけ

+ vertex.color=ifelse(V(g)$Faction==1,"Pink","Lightgreen"),

+ vertex.label.color="gray50", #ノードのラベルの色

#ノードのラベルのスタイル 1: 普通, 2: 太字, 3: 斜体, 4: 太字斜体, 5: ギリシャ文字

+ vertex.label.font=2,

+ vertex.frame.color="white", #ノードの枠の色

+ vertex.label.cex=0.8, #ノードラベルの文字サイズ

+ edge.width=E(g)$weight, #エッジ属性weightをエッジの太さとする

+ edge.color="gray80", #エッジの色

+ layout=layout.fruchterman.reingold) #ネットワークのレイアウト手法

RColorBrewerを使ってノードをカテゴリーごとに美しく色分け

ノードをカテゴリー（例えばGene Ontologyなど）ごとに色分けする場合、ColorBrewerを利用したい場合があると思います。少し面倒ですが、色とカテゴリーの対応表を作成することで可能になります。具体的には以下のようです。ここではノード10個で構成されるネットワークを例に考えます。

> library(RColorBrewer) #ColorBrewerの読み込み

> g <- graph.ring(10) #ネットワーク

# ノードのカテゴリー（一般にはファイルから読み込む）

> class <- c("Energy metabolism","Energy metabolism","Energy metabolism",

+ "Lipid metabolism","Lipid metabolism","Lipid metabolism",

+ "Lipid metabolism","Cell cycle","Cell cycle","Cell cycle")

# ColorBrewerのパレットを設定して（ここではBrBG）

# カテゴリーと色の対応表を作る。

> clscolor <- data.frame(class=levels(as.factor(class)),

+ color=I(brewer.pal(nlevels(as.factor(class)),"BrBG")))

#ネットワークの描画

> plot(g,

+ vertex.color=clscolor$color[match(class,clscolor$class)],

+ vertex.label=class)

値で色分けする場合：Gene Ontologyのような機能クラスではなく、遺伝子発現量のような連続値でノードを色分けしたい（例えば、値が高いものを濃い色、値が低いものを薄い色で描画する）場合はその値の範囲ごとにクラス分けをすれば同様に色付けすることができます。例は以下のようです。

> g <- graph.ring(10) #ネットワーク

# ノードの値（実際にはファイルから読み込み）

> val <- rnorm(10)

# 数値をクラス分け（ここでは8分割）

# 使用するパレットに応じて変更

> class <- cut(val,breaks=8)

# ColorBrewerのパレットを設定して（ここではYlOrRd）

# カテゴリーと色の対応表を作る。

> clscolor <- data.frame(class=levels(as.factor(class)),

+ color=I(brewer.pal(nlevels(as.factor(class)),"YlOrRd")))

# ネットワークの描画

> plot(g,

+ vertex.color=clscolor$color[match(class,clscolor$class)],

+ vertex.label=val)

レイアウトアルゴリズムを使ってネットワーク描画

Kamada-Kawaiのアルゴリズムの場合

plot(g,layout=layout.kamada.kawai)

Fruchterman-Reingoldのアルゴリズムの場合

plot(g,layout=layout.fruchterman.reingold)

自動（最適なアルゴリズムを勝手に適応して）で描画

plot(g,layout=layout.auto)

ネットワークの描画に関しては Static and dynamic network visualization with R がとても参考になります。静的描画はもちろんのこと動的描画も充実しています。

各ノードのネットワーク特徴量を出力

次数（近傍ノードの数）

> degree(g)

クラスタ係数

良く知られた Watts & Strogatz (1998) の平均クラスタ係数は以下の様にして求めることができます。

> transitivity(g,type="local",isolates="zero")

Barrat et al. (2004) による重み付きネットワークにおけるクラスタ係数は以下の様に求められます。

> transitivity(g,vids=V(g),type="weighted",isolates="zero")

Centrality 中心性 (Freeman, 1979)

「どのノードが重要か」はネットワーク科学において重要です。特に、Centralityは様々な特徴と相関があり、ノードの性質を理解、予測する上で役に立ちます。基本的に、重みやエッジの方向性はattributeとして与えられていれば自動的に考慮されます。

次数中心性（次数）

> degree(g) / (vcount(g) - 1)

媒介中心性 (Betweenness centrality)

> betweenness(g)

近接中心性 (Closeness centrality)

> closeness(g)

固有ベクトル中心性 (eigenvector centrality)

> evcent(g)$vector

Google Page Rank

Page rankの計算に使われるアルゴリズムのひとつ。これだけで、Google Page Rankを計算している訳ではないことに注意。

> page.rank(g)$vector

重み付きの場合に関しては g に重み情報が割り振られている場合（E(g)$weight）、デフォルト重み付きで計算されます。

ネットワーク全体の特徴量

次数分布P(k)

> degree.distribution(g)

> plot(degree.distribution(g)) # プロットする場合。

※ R ver. 3.0.0以上 + igraph ver. 0.6.5以上である場合、degree.distribution関数は利用不可能であるようです。取りあえず、以下の様にすると次数分布を得ることができます。R ver. 3.0.3 + igraph ver. 0.7.0で確認したところ問題は解決されていました。以下の様にしても、次数分布を得ることができます。

> table(degree(g))/sum(table(degree(g)))

> plot(names(table(degree(g))),table(degree(g))/sum(table(degree(g)))) #プロットする場合

P(k)を応用すると次数エントロピー（-\sum_{k=1}^{N-1} P(k) log P(k)）も計算できます。

> dist <- table(degree(g))/sum(table(degree(g)))

> -sum(dist*log(dist))

P(k)がベキ分布（i.e., k^{-r}）であると仮定した場合の指数 r の推定法 (Clauset et al., 2007)

> 1+vcount(g)/sum(log(degree(g)/min(degree(g))))

クラスタ係数

Newman et al., (2002) による定義

> transitivity(g,type="global")

Assortative Coefficient

連結する二ノード間の次数の相関係数（Newman, 2002）

> assortativity.degree(g)

平均最短経路長

> average.path.length(g)

各ノード間の最短経路長を行列で表示

> shortest.paths(g)

ノードAからBの最短経路を全て列挙する。

> get.all.shortest.paths(g,"A","B",mode="out")

最短経路長の頻度

> path.length.hist(g)

ネットワークの処理

最大 (弱) 連結成分の抽出

> cls <- clusters(g,"weak")

> g_comp <- delete.vertices(g,subset(V(g),

+ cls$membership!=which(cls$csize==max(cls$csize))[[1]]))

ヌル仮説ネットワーク（ネットワークのランダム化）

ネットワークの構造指標が有意であるかどうかを検証する際には、次数を固定したままランダム化されたネットワークとの比較を行う場合があります（例えば、[Milo et al, 2002; 2004]参照）。

単純なエッジ交換アルゴリズム（次数を保存したランダム化ネットワークを作成する方法）を用いる場合はrewire関数を用います。多重エッジや自己ループができないようにkeeping_degseqのオプションを与えておきます。ただ、この関数は論文（Milo et al., 2002; 2004）で用いられたアルゴリズムとは異なり、相互エッジ（<=>）と単方向エッジ（=>）の数を保存しないので、注意が必要です。niterはエッジ交換の回数を意味します。全てのエッジが交換されることを終了条件と考えると、これはクーポンコレクター問題と（ほぼ）同じになります。そのため、例えば、エッジ数を E とすると、E x (log E + c) / 2 回エッジ交換を行った場合（2で割られているのは1回のエッジ交換で2本のエッジが交換されるため）、各エッジが少なくとも1回は交換されている確率は近似的に exp(-exp(-c)) となります。ここで c は定数です。

> g_rand <- rewire(g, with=keeping_degseq(niter=ecount(g)*(50+log(ecount(g)))/2,loops=F))

似た関数にrewire.edgeがありますが、これは均一に選択されたエッジが均一に選択されたノードに張り替える関数なので、ランダム化には不適切である（ノード次数を保存しない）ので注意。

他にも、Configuration model（例えば、（Catanzaro et al., 2005）参照）から得られるランダムネットワークをヌル仮説ネットワークとして用いる場合もあります。

> deg <- degree(g) #次数情報

#次数に基づいてランダムネットワークを作成（無向グラフの場合）

> g_rand <- degree.sequence.game(deg)

#連結性を保証し、多重エッジの生成を回避する場合（推奨）

> g_rand <- degree.sequence.game(deg,m="vl")

#有向グラフの場合

> deg_out <- degree(g,m="out") #出次数

> deg_in <- degree(g,m="in") #出次数

> g_rand <- degree.sequence.game(deg_out,deg_in,m="vl") # 推奨

# ただし、条件によって実行できない場合がある。その場合は m="vl"を除く

デフォルト状態において、この関数から生成されるネットワークは連結性が保証されず、多重エッジの存在が許されます。オプションm="vl"をつけることで、ランダム化されたネットワークの連結性が保証され、エッジの交換で多重エッジの生成を避けることができます。ただし、元になるネットワークの連結性が保証されていないといけないことに注意してください。例えば、最大連結成分を入力にすればOKです。ただ、この関数も、相互エッジと単方向エッジの数を保存しないという点で問題になることがあるので使う際には注意が必要です。また、vlモードでネットワークを作成できない場合があります。この際、オプションm="vl"を除くことで実行できるようになりますが、多重エッジや自己ループの生成を許すので注意が必要です。

相互エッジの数、単方向エッジの数、入次数、出次数を保存し、多重エッジや自己ループが生成されないような（つまり元論文と同じ）アルゴリズムを使いたい場合は、edge_swapping_randomization関数を用いてください。

> source("edge_swapping_randomization.R")

> g_rand <- edge_swapping_randomization(g)

無向、有向ネットワークのどちらでもつかえます。ただ、重み付きネットワークには対応していません。

重み付きの無向ネットワークをランダム化したい場合は、edge_swapping_randomization_weighted_undirected関数を用いてください。エッジを交換する際、重みも一緒に交換されるようなランダム化です。

コミュニティ検出 (Non-overlapping)

ネットワークをグループ分け（コミュニティ分割）することはいろいろな局面で重要です（Fortunato, 2010）。ここでは、Non-overlapping community detection（ノードがひとつのみのコミュニティに属する場合のコミュニティ抽出）を行います。

以下では、結果をdataに格納することを考えます。

Edge betweennessに基づくアルゴリズム (Girvan-Newman algorithm)

（Girvan & Newman, 2002）参照

data <- edge.betweenness.community(g)

ランダムウォークに基づくアルゴリズム

（Pons & Latapy, 2006）参照

> data <- walktrap.community(g, modularity=TRUE)

貪欲アルゴリズムに基づく高速なアルゴリズム

（Clauset et al., 2004）参照

> data <- fastgreedy.community(g)

固有ベクトルに基づくアルゴリズム

（Newman, 2006）参照

> data <- leading.eigenvector.community(g,options=list(maxiter=1000000, ncv=5))

多段階最適アルゴリズムに基づく

（Blondel et al., 2008）参照

> data <- multilevel.community(g)

スピングラスに基づくアルゴリズム

（Reichardt & Bornholdt, 2006）参照

> data <- spinglass.community(g)

全探索（全てのコミュニティアサインメントを試すこと）による方法

計算時間がかかるため、小さなネットワークのみ有効

> data <- optimal.community(g)

ラベル伝播法に基づくアルゴリズム

（Raghavan et al., 2007）参照

> data <- label.propagation.community(g)

Infomap法に基づくアルゴリズム

（Rosvall & Bergstorm, 2008）参照

> data <- infomap.community(g)

各ノードのメンバーシップ（どのノードがどのコミュニティに属しているかの情報）の参照

> data$membership

最適分割時、のQ値（モジュラリティー）とコミュニティの数

> max(data$modularity)

> max(data$membership)

任意の計算ステップにおけるコミュニティ分割のメンバーシップを得るためには以下の関数を使います。関数dendPlotを使うと、ステップ数とクラスタリングの階層関係が分かりますので、この情報から適切なステップ数（steps=のパラメータ）を選択してください。なお、この関数が使えるのは、階層的にクラスタリングする以下のアルゴリズムのみで有効です。

• edge.betweenness.community

• walktrap.community

• fastgreedy.community

> community.to.membership(g,data$merges, steps=2)

どの手法を選択すればよいのでしょうか？：計算時間と検出精度のトレード・オフ

さて、様々な手法が利用できるのは良いことなのですが、結局のところ、どれを使えばよいのでしょうか。それぞれの手法には利点と欠点があります。目的と状況に応じて選択することが重要です。

基準になるのは計算時間（とにかく早く計算したい）と検出精度（より最適な解をを得たい）だと思います。Fortunato (2010) はコミュニティ検出の多くの手法について、計算時間と検出精度の視点から比較しています。

計算時間では、貪欲法（fastgreedy.community）が最短です。計算時間が短い順に挙げると次のようになります：貪欲法 < Girvan-Newman法 < 固有ベクトル法 < 焼きなまし法。

検出精度（どれだけ適切にモジュールを見つけられているか）では、焼きなまし法（spinglass.community）が最高です。精度の高い順に挙げると次のようになります：焼きなまし法 > 固有ベクトル法 > 貪欲法 > Girvan-Newman法。ただ、この結果は、あるベンチマークを用いた検証から得られたものですので、解釈には注意が必要でしょう。あくまでひとつの基準として考えて下さい。

この結果から分かるように、計算時間と検出精度はトレード・オフの関係にあります。これは手法を選択するひとつの判断基準になるでしょう。

コミュニティ抽出 (Overlapping)

ノードが複数のコミュニティに属することを許した場合のコミュニティ抽出を行います。

これはlinkcommというパッケージを用います（igraphではないことに注意してください。その他のパッケージも一部読み込ますが、読み込みは自動的に行われるので気にする必要はありません）。

> library(linkcomm)

ノード名が文字列で与えられてるエッジリストを入力とします。例えば以下のようです。最後の数字はエッジの重みに対応します。

A B 1

B C 2

C D 3

...

ネットワークデータとして読み込みます。

> d <- read.table(filename)

> gw <- d # そのまま重み付きネットワークとして認識される。

> g <- d[1:2] # 2行目までをつかえば重みなしネットワーク

Link community（Ahn et al., 2010）の場合

無向ネットワークの場合（重みが付いている場合（gw）は自動的に重みが考慮される）

> data <- getLinkCommunities(g)

有向ネットワークの場合（重みが付いている場合（gw）は自動的に重みが考慮される）

> data <- getLinkCommunities(g,directed=T)

コミュニティ抽出の結果をネットワークで表示する。

> plot(data,type="graph")

エッジ1本のみで構成されるコミュニティは表示されないことに注意して下さい。

typeの与え方によって様々な表示が可能

summary: 分割密度（D値）とデンドログラムを表示

members: コミュニティのメンバーシップ行列を表示

commsumm: 各コミュニティのモジュラリティを棒グラフで表示

dend: デンドログラムを各コミュニュティごとに色付きで表示

任意の閾値（D値）でコミュニティ抽出を行う。

> lc <- newLinkCommsAt(data,cutat=0.95)

ノードのメンバーシップ情報を出力

> lc$nodeclusters

エッジ1本のみで構成されるコミュニティの情報は出力されないことに注意して下さい。

Overlapping Cluster Generator (OCG)（Backer et al., 2012）の場合

使い方は基本的にgetLinkCommunitiesと同じです。以下のように行います。

> oc <- getOCG.clusters(g)

plotはgraphとmembersが利用可能。

> plot(oc,type="graph")

ネットワークモチーフを見つける

Configuration modelをヌル仮説ネットワークとして、それと比較して有意に頻出するサブグラフを見つけます。ただし、これはNetwork motifの元論文（Milo et al., 2002）で用いられた手法とは異なることに注意してください（ネットワークのランダム化の手法が異なります。ネットワークのランダム化も参照）。

実際のネットワークのサブフラフを数える。以下は3-node subgraphの場合を考える。4-node subgraphの場合は3を4に置き換えるだけでよい。ただし5以上は対応していない。

c_real <- graph.motifs(g,3)

ヌル仮説ネットワーク（無向）を1000個作成して、それぞれのサブフラフの数をc_randに格納する。

deg <- degree(g)

c_rand <- data.frame()

for(i in 1:1000){

g_rand <- degree.sequence.game(deg)

c_rand <- rbind(c_rand,graph.motifs(g_rand))

}

names(c_rand)<-0:3

有向ネットワークの場合は以下のようにする。

deg_out <- degree(g,m="out")

deg_in <- degree(g,m="in")

c_rand <- data.frame()

for(i in 1:1000){

g_rand <- degree.sequence.game(deg_out,deg_in)

c_rand <- rbind(c_rand,graph.motifs(g_rand,3))

}

names(c_rand)<-0:15

各サブグラフのZ-scoreを出力

例）

無向ネットワークの場合

3: Triad (Undirected case)

有向ネットワークの場合

7: Feedforwaed loop, 11: Feedback loop

> (c_real-colMeans(c_rand))/sqrt(colMeans(c_rand**2)-colMeans(c_rand)**2)

モチーフの種類とIDの対応は以下のようにして確認できる。

> plot(graph.isocreate(3,7))

任意のモチーフのIDを得る場合は以下のようにして確認できる。

任意のグラフを入力してg_sに確認したとする。

> graph.isoclass(g_s)

ネットワークに埋め込まれたモチーフを構成するノード番号を出力する。 以下はFFLの場合（幾何的な制約を考慮していないので、3-node motifの場合、graph.motifsでカウントされた数より大きくなる）。

> graph.get.subisomorphisms.vf2(g,graph.isocreate(3,7))

ネットワーク比較

和

> graph.union(g1,g2)

積

> graph.intersection(g1,g2)

差分

> graph.difference(g1,g2)

ただし、|g1| > |g2|

Box-covering algorithm：ネットワークのフラクタル次元

現実のネットワークはフラクタル（自己相似）構造を示す場合がしばしばです（Song et al., 2005）。このようなフラクタル構造はネットワーク進化の理解や機能モジュールの同定において重要な洞察を与えてくれます。

ネットワークのフラクタル性は、Box-covering algorithmによって判定されます。具体的に、様々な解像度（l_B）でネットワークを粗視化した際、その解像度と被覆に用いたBoxの数の統計的な傾向から特徴付けられます。これは、ボックスカウント法によるフラクタル性の判定法のネットワーク版だと捉えることができます。

Box coveringはグラフ彩色問題と関連があり、厳密解を得る場合膨大な計算時間を必要とします。そのため、近似アルゴリズムが役立ちます。Song et al. (2007) はさまざまな近似アルゴリズムを紹介しています。

最も実装が簡単なのは貪欲法に基づくアルゴリズムです。あるネットワークGに対する解像度（l_B）におけるBox coveringは、G から変換されたネットワーク G'（このネットワークにおいて頂点 i と j は G における最短距離が l_B 以上であれば隣接している）の彩色問題として捉えることができます。これは、貪欲法によるグラフ彩色を応用して解くことができます。

以下に、貪欲法によるBox coveringのRのスクリプト（Greedy Box-Covering Algorithm）の例を示します。最適解（最少数の色でのグラフ彩色）を効率的に得るためには、ノードの順序（idの割り振り方）などに工夫の余地が残されています。以下を参考に改良してもよいでしょう。

# ネットワークオブジェクト

g <- barabasi.game(20,1,1,directed=F)

l_B_max <- diameter(g) + 1

node_num <- vcount(g)

pathlength <- shortest.paths(g)

color <- matrix(1,ncol=l_B_max,nrow=node_num)

node_id <- 1:node_num

# 貪欲法によるBox covering

for(i in 2:node_num){

for(l_B in 1:l_B_max){

data <- as.data.frame(cbind(node_id,pathlength[i,]))

neighbor <- subset(data,data[[1]] < i & data[[2]] >= l_B)

unused_col <- setdiff(color[,l_B],color[,l_B][neighbor[[1]]])

if(length(unused_col) > 0){

idx <- as.integer(runif(1,max=length(unused_col))+1)

color[i,l_B] <- unused_col[idx]

}

else{

color[i,l_B] <- max(color[,l_B]) + 1

}

}

}

# 例えばl_B=2の場合のメンバーシップ（被覆状態）を見る場合は以下のよう。

color[,2]

# フラクタル性の判定（単純に、ベキ関数と指数関数のどちらでよくフィットできるかで判定）

# データの準備

data <- data.frame()

for(l_B in 1:l_B_max){

data <- rbind(data,c(l_B,max(color[,l_B])))

}

names(data) <- c("l_B","N_l_B")

# ベキ関数フィッティング

sum_pow <- lm(log(N_l_B)~log(l_B),data=data)

# 指数関数フィッティング

sum_exp <- lm(log(N_l_B)~l_B,data=data)

# 判定

if(summary(sum_exp)$fstatistic[[1]] > summary(sum_pow)$fstatistic[[1]]){

cat("This network is NOT fractal.\n")

}

else{

cat("This network is fractal.\n")

# フラクタル次元の表示

cat("Fractal dimension d_B: ",-sum_pow$coefficients[[2]],"\n")

}

Song et al. (2007) も言及していますが、この貪欲法によるBox covering algorithmによって被覆されたBoxの連結性は保証されていません。例えば、2種類のBoxを考えると（○-●-○）のような分割が出現する場合があります。Boxの連結性を保証したい場合はMaximum-excluded-mass-burning (MEMB) algorithmを用いる必要があります。

Bow-Tie Decomposition: Finding Input and Output Layers of Networks

現実のネットワーク（辺に方向性有り）はしばしばちょうネクタイ（Bow tie）構造を有します。具体的に、そのようなネットワークには「入力層・コア・出力層」が存在します（辺に向きが無い場合は、コア・周辺構造（Core-Periphery Structure）と呼ばれることもあります（Csermely, 2012））。

特に生物学において、このBow-tie構造はロバストネスと関連があり、環境適応、ノイズ軽減、Evolvabilityなどの文脈で語られ（Kitano, 2004）、生命システムの理解に役立ちます。

Bow-Tie Decompositionの単純な手法としては強連結成分分解に基づくものがあります（Yang, 2011）。ここでは、最大の強連結成分がコアとして定義されます。入力層と出力層は、それぞれ、コアに到達できるがコアからは到達できないノード集合、コアからは到達できるがコアには到達できないノード集合として定義されます。

以下に、有向ネットワークからコア・出力層・入力層を見つけだすスクリプトの例を挙げます。チューブ（入力層と出力層を直接つなげる経路）やひげ（入力／出力層からコア以外にのびる経路）なども、これを応用して、見つけ出すことができます。

library(igraph) #igraphパッケージの読み込み

g <- erdos.renyi.game(20,0.08,directed=T) #グラフオブジェクト

# もし、ノードラベルが「与えられていない」なら

if (length(V(g)$name) == 0) {

V(g)$name <- 1:vcount(g)

}

cls <- clusters(g,"strong") #強連結成分分解

#最大強連結成分（LCC）に属するノード

lcc_nodes

<- subset(V(g)$name,cls$membership==which(cls$csize==max(cls$csize)))

#入力層・出力層の候補（LCCに属さないノード）

candidate

<- subset(V(g)$name,cls$membership!=which(cls$csize==max(cls$csize)))

#各ノードの最短距離を計算

dist <- shortest.paths(g,m="out")

dist_from_lcc <- dist[lcc_nodes,]

node_name_from_lcc <- colnames(dist_from_lcc)

#LCCから到達できないノード集合

unreachable_from_lcc

<- node_name_from_lcc[unique(which(dist_from_lcc==Inf,arr.ind=T)[,2])]

#LCCから到達できるノード集合

reachable_from_lcc

<- node_name_from_lcc[unique(which(dist_from_lcc!=Inf,arr.ind=T)[,2])]

#入力層と出力層のノードを探す

substrate <- c() #入力層ノード集合

product <- c() #出力層ノード集合

for (cand in candidate){

#LCCに到達できるが、LCCからは到達できないノードが入力層ノード

reachable_from_cand <- names(subset(dist[cand,],dist[cand,]!=Inf))

if (length(intersect(reachable_from_cand,lcc_nodes)) > 0

&& is.element(cand,unreachable_from_lcc) == T) {

substrate <- c(substrate,cand)

}

#LCCから到達できるが、LCCには到達できないノードが出力層ノード

unreachable_from_cand <- names(subset(dist[cand,],dist[cand,]==Inf))

if (length(intersect(unreachable_from_cand,lcc_nodes)) > 0

&& is.element(cand,reachable_from_lcc) == T) {

product <- c(product,cand)

}

}

#結果の出力

#入力層・出力層・コアのメンバーシップを作る。

mem_sub <- ifelse(is.element(V(g)$name,substrate)==T,1,0)

mem_prod <- ifelse(is.element(V(g)$name,product)==T,2,0)

mem_lcc <- ifelse(is.element(V(g)$name,lcc_nodes)==T,3,0)

mem <- mem_prod + mem_lcc + mem_sub

#ネットワークの描画

plot(g,vertex.color=mem)

Threshold selection based on random matrix theory for constructing correlation networks：ランダム行列理論に基づく閾値選択で相関ネットワークを構築

遺伝子の共発現、（微）生物の共起関係、脳部位の機能関係を考える場合、相関行列からネットワーク（相関ネットワーク）を得たい場合があります。一般的に、相関係数の閾値を用いて二値化することで、ネットワークを得ます。ただ、どのように閾値を設定すれば良いのでしょうか。統計指標に基づく方法やネットワーク指標に基づく方法もありますが、ここでは、最近比較的よく使われている、ランダム行列理論（RMT）に基づく閾値選択について紹介します。

RMTに基づく手法では、閾値処理した相関行列の固有値「間隔」の分布に基づいて最適な閾値を選択します。具体的に、固有値がなるべく相関しない（固有値「間隔」の分布が指数分布に最も近くなる）ような行列を得ることを考えます（なお、固有値が相関する場合は固有値の間隔はWigner–Dyson分布に従います）。これは（誤解を恐れずに、単純化して言うならば）ある行列からランダム性を排除し、システム特異な（非ランダムな）性質を取り出す作業に対応します。

以下は、ランダム行列理論に基づく閾値選択のための関数 thresholding.RMT のソースコードです（もちろん igraph のライブラリを要求します）。

thresholding.RMT <- function(mtx){

mtx <- abs(mtx)

rc_opt <- 0

ks_opt <- 1

for(rc in seq(0.300,0.999,by=0.001)){

mtx_tmp <- ifelse(mtx > rc,mtx,0)

ei <- eigen(mtx_tmp)

diff_eigenvalue <- diff(sort(ei$values))

mean_diff_eigenvalue <- mean(diff_eigenvalue)

if(mean_diff_eigenvalue > 0){

diff_eigenvalue <- diff_eigenvalue / mean_diff_eigenvalue

# calc. Kolmogorov-Smirnov distance

res.ks <- ks.test(diff_eigenvalue,"pexp")

if(res.ks$statistic < ks_opt){

rc_opt <- rc

ks_opt <- res.ks$statistic

p.value_opt <- res.ks$p.value

}

}

}

cat("Optimal threshold value = ",rc_opt,"\nKS distance = ",ks_opt," (p = ",p.value_opt,")\n",sep="")

mtx_opt <- ifelse(mtx > rc_opt,1,0)

g <- graph.adjacency(mtx_opt,mode="undirected",weighted=NULL)

return(g)

}

相関行列を与えれば、最適な閾値から作成されるネットワーク（グラフオブジェクト）を返します。なお、KS距離が小さい（P値が大きい）ほど指数分布への当てはまりが良いことを意味します。

# 仮想的な（100x100の）相関行列（対称行列）の作成

> mtx <- matrix(runif(100*100),ncol=100,nrow=100)

> mtx <- (mtx + t(mtx)) / 2

> diag(mtx) <- 1 # 自分自身との相関係数は必ず1になるので

# 最適な閾値の計算（出力結果は例です）

> g <- thresholding.RMT(mtx)

Optimal threshold value = 0.884

KS distance = 0.08415366 (p = 0.4847843)

Finding Driver Nodes (Minimum Dominating Set) Using Linear Programming

線形計画法を用いて最小支配集合を求めるスクリプトです。RパッケージlpSolveが必要です。

あるシステムを反映するネットワークにおいて、そのシステムのダイナミクスを支配するノード（Driver node）を見つけだすことは、制御の文脈において極めて重要です。

そのようなDriver nodeを見つけるための最も素朴なアプローチは、最小支配集合を求めることでしょう。事実、最小支配集合は生物学や工学分野で幅広く用いられています（Nacher and Akutsu, 2012）（ただ、これは Lui et al. (2011) が示すDriver nodeとは厳密には異なるので注意）。

ところが、最小支配集合を正確に求めることは計算量の側面からとても難しい問題であると考えられています。このような場合、最適化手法が有効です。

ここでは線型計画法を用いて、最小支配集合を求めてみたいと思います。

ノード i が最小支配集合に属するか否かは、変数 x_i で記述されるとします。つまり x_i = 1 なら最小支配集合に属する、x_i = 0なら属さないことを意味します。

このとき、関数 f(x) を最小化することを考えます。N個のノードで構成されるネットワークの場合以下のようになります。

f(x) = x₁ + ... + x_N

ただし、以下のような制約を考えます。具体的に、すべてのノード（i = 1, ... , N）に対して以下が成り立つように関数 f(x) を最小化させます。

x_i + ∑_j=1^N A_ij x_j ≥ 1

ここで A_ij はネットワークの隣接行列で、j => i という辺が存在するなら 1、そうでなければ 0 となります。

具体的に、以下のスクリプトで最小支配集合を求めることができます。

#igraphパッケージの読み込み

library(igraph)

#線型計画法LPパッケージの読み込み

library(lpSolve)

#グラフオブジェクト

g <- erdos.renyi.game(20,0.08,directed=T)

#LPのための下準備

m_g <- t(as.matrix(get.adjacency(g)))

diag(m_g) <- 1

f.dir <- rep(">=",vcount(g))

f.rhs <- numeric(vcount(g)) + 1

f.obj <- numeric(vcount(g)) + 1

#0-1整数計画法で最小支配集合を求める（緩和したいならall.bin=Tを除く）

res <- lp("min", f.obj, m_g, f.dir, f.rhs, all.bin=T)

#ノードのメンバーシップ

res$solution

# 1: 最小支配集合に属する

# 0: 属さない

#ネットワーク描画

plot(g,vertex.color=res$solution)

最小支配エッジ集合は、Line graphの最小支配頂点集合から求めることができます。Line graphは関数 line.graph で求めることができるので、応用として挑戦してみてください。

Finding Maximum Matching Using Linear Programming

線形計画法を用いて、ネットワークの最大マッチングを求めます。lpSolveのパッケージが必要です。

最大マッチングはその名の通り、あらゆる場合における「マッチング」を考える際に役立ちます。2012年のノーベル経済学賞はこのマッチングに関する研究に送られたので、ご存知の方も多いかと思います（例えばここが参考になります）。

また、システムの制御可能性（Controllability）の側面から言及される場合（Liu et al., 2011）もあり、ネットワークのマッチングは幅広い応用が期待できます。

ここでは、以下のような整数計画問題におきかえて、最大マッチングを求めることを考えます。

頂点数 N と辺数 E で構成されるネットワークにおいて、辺 i が最大マッチングに属しているか否かは変数 x_iで表します：x_i = 1（辺 i が最大マッチングに属する）、x_i = 0（属さない）。

この時、以下の関数を最大化します。

f(x) = w₁x₁ + ... + w_Ex_E

ここで、w₁, ... , w_E は辺の重みです。

制約条件としては、以下を考えます。マッチングは、全てのノード（i = 1, ... , N）に対して、マッチングに属する辺が高々 1 本しかないわけですから、以下が成り立つようにします。

∑_j=1^E B_ij x_j ≤ 1

ここで、B_ij はノード i が辺 j の端点であるなら 1、そうでないなら 0 を示す行列です。

具体的に、以下のスクリプトを用いて、最大マッチングを求めることができます。

#igraphパッケージの読み込み

library(igraph)

#線型計画法パッケージの読み込み

library(lpSolve)

#グラフオブジェクト

g <- erdos.renyi.game(20,0.2,directed=F)

#辺の重みを与える。

E(g)$weight <- numeric(ecount(g)) + 1

# ノードの名前を与える（すでに与えられている場合は必要ない）

V(g)$name <- 1:vcount(g)

#LPのための下準備：B_ijの作成

edge_list <- get.edgelist(g)

n_col <- nrow(edge_list)

f.con <- matrix(0,nrow=vcount(g),ncol=n_col)

for(i in 1:vcount(g)){

idx <- which(edge_list[,1] == V(g)$name[[i]] | edge_list[,2] == V(g)$name[[i]])

f.con[i,idx] <- 1

}

f.dir <- rep("<=",vcount(g))

f.rhs <- numeric(vcount(g)) + 1

f.obj <- E(g)$weight

#0-1整数計画法で最大マッチングを求める（緩和したいならall.bin=Tを除く）

res <- lp("max", f.obj, f.con, f.dir, f.rhs, all.bin=T)

#エッジのメンバーシップ

res$solution

# 1: マッチング集合に属する

# 0: 属さない

#ネットワーク描画

plot(g,edge.color=res$solution+1,vertex.size=5,edge.width=5)

Functional Cartography of Complex Networks

大規模なネットワークにおいてノードの役割や性質が分類できると便利ですよね。ひとつの考え方としてFunctional Cartography（Guimera & Amaral, 2005）という分類法があります。これはコミュニティ検出アルゴリズムなどをもちいてネットワークをモジュール分割し、モジュール内・外に対してノードがどのような役割を果たしているのかを特徴付けます。具体的には以下のように分類されます。

• R1: Ultra-peripheral nodes. 「超」周辺ノード

• R2: Peripheral nodes. 周辺ノード

• R3: Non-hub connectors. ハブではないが、複数のモジュールをつなぐコネクター

• R4: Non-hub kinless nodes. ハブではないが、多くのモジュールにつながるノード

• R5: Provincial hubs. モジュール内のハブ

• R6: Connector hubs. 複数のモジュールをつなぐハブ

• R7: Kinless hubs. 多くのモジュールにつながるハブ

例えば、ハブの存在するネットワークではハブの攻撃に対して脆弱（Albert et al., 2000）ですが、モジュール内のハブを攻撃するのと複数のモジュールをつなぐハブを攻撃するのとでは挙動が違います（Han et al., 2004）。このように、この分類で重要そうなノードを絞り込むことができます。

計算するRスクリプトは、ちょっと長くなったので、ここからダウンロードして下さい。The R script is downloadable here.

なお、Functional cartographyはrnetcartoパッケージでも計算可能です。こちらの方を使った方が簡単でしょう。