【R】ネットワーク解析を学んでみたい④

・お題：生体内の様々な分子は相互作用し、絶妙なバランスを保っている。生体内の分子の役割を理解するのに、ネットワーク解析が有効らしい。学んでみたい。

・私は素人なので、参考サイトをなぞるような形で学んでいきたい。今回は以下のサイトを参考に学ばせていただく。九州工大の竹本先生のHP。日本語でここまで詳しく解説されているHPは他に見たことがないので、ぜひこちらのHPをご覧いただきたい。

sites.google.com

・竹本先生の貴著「生物ネットワーク解析」も名著だと思うので、気になる方はぜひチェックしてみてほしい。

www.coronasha.co.jp

・前回、コミュニティ抽出をやってみた。今回はネットワークモチーフを見つけてみたい。ネットワークモチーフは、あるネットワークの中で多く観測される部分ネットワークのことで、何等か機能的意義のある構成要素である場合があるらしい。多い少ないの判断は統計的な優位性で判断するらしい。統計……かぁ…（遠い目）。なにはともあれ、やってみる。

・まず、モチーフを調べたいグラフを準備する。以下をtextファイルに打ち込んでカレントディレクトリに保存した。

A B 1
B C 1
C A 1
C D 1
D A 1
D E 1

・読み込んで図示。

> library(igraph)
> library(RColorBrewer)
> d <- read.table("test.txt")
> g <- graph.data.frame(d[1:2],directed=F)
> plot(g)

・あらかじめ設定されたモチーフが何個あるか数え上げて、それがランダムに設定されたネットワークと比べて多く出現するか否か調べる感じになる。

・「あらかじめ設定されたモチーフ」が何なのかというと、例えば3-node subgraphがあるか見るのであれば、3つのノードでできているグラフのパターンを考えてやれば良いみたい。無向グラフなら、パターンは三角形の辺の数（0,1,2,3本）と同じく、4種類しかない。

・求めてみると、以下のような感じ。エッジの無いノードがないので、最初と2番目のモチーフはそもそも想定されないみたい。3番目のモチーフは4つ、4番目のモチーフは2つ見つかった。
> graph.motifs(g,3)
[1] NA NA 4 2

・ちなみに、4つのノードでできるグラフのパターンは11個あるらしい。求めてみるとこんな感じ。

> graph.motifs(g,4)
[1] NA NA NA NA 0 NA 2 1 0 1 0

・どれがどんなモチーフに対応しているかはgraph.isocreate関数で描くことができる。

> plot(graph.isocreate(4,7,directed=F))

・次に、もう少し大きいネットワークを考える。以下の文字列をtextで保存し、先ほどと同様に読み込み（割愛）、図示する。

A C 1
D C 1
C B 1
E C 1
C G 1
C H 1
G F 1
G I 1
H I 1
I Q 1
Q R 1
R M 1
R K 1
J K 1
K L 1
K M 1
L M 1
M O 1
M P 1
G H 1
F I 1
O R 1

>plot(g)

・これに対して、3つのノードのグラフモチーフを調べてみる。
> c_real <- graph.motifs(g,3)
> c_real
[1] NA NA 38 6

・次に、このグラフの次数を算出する。

> deg <- degree(g)
> deg
A D C E G H I Q R J K L M F O B P
1 1 6 1 4 3 4 2 4 1 4 2 5 2 2 1 1

・この次数を保存したまま、ランダムなグラフを描いてみる。これは実行するたびにグラフが変わる。

> plot(degree.sequence.game(deg))

・このランダムなネットワークを1000個作成し、それぞれのグラフに対して先ほどの3ノードモチーフがどんな感じに含まれていたかを調べる。

> 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

> c_rand
0 1 2 3
1 NA NA 20 1
2 NA NA 41 2
3 NA NA 40 0
4 NA NA 36 0
5 NA NA 35 2
6 NA NA 36 2
7 NA NA 34 3
8 NA NA 44 4

…（1000行続くことになるが、表示は端折られる）

・ランダムに作成したネットワークに対して、自前で用意したネットワークの方が、モチーフをたくさん含んでいたりすると、モチーフに意義があるような気がする。。ここでたくさんか否かを判断する基準として、統計を用いるっぽい。Z-scoreを求めてみる。

・Z-scoreは、気になる値から平均値を引いて標準偏差で割ることで求められる。今回でいうと、気になる値はc_real、平均値はcolMeans(c_rand)になる。標準偏差は分散の平方根なので、sqrt(colMeans(c_rand**2)-colMeans(c_rand)**2)になる。

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

0 1 2 3
NA NA 0.6080921 3.6046358

・ちなみに、c_randの2番目のモチーフは以下のような分布。私の用意したグラフgは、有意にこのモチーフを多く含んでいたわけではなさそう。。

>hist(c_rand$`2`)

・最後に、グラフの重ね合わせをしていったん終わりたいと思う。その先は今の私には厳しいので、ぼちぼち勉強していきたい。。

・以下のようにtest1とtest2のテキストファイルを作成し、カレントディレクトリに保存した。

test1.txt

A B 1
B C 1
C D 1
D A 1
E A 1
F E 1
G E 1
H E 1
I E 1

test2.txt

A B 1
B C 1
C D 1
D A 1
B D 1
J D 1
K D 1
J K 1
K L 1
K M 1

・それぞれを読み込んで図示してみる。

> d1 <- read.table("test1.txt")
> g1 <- graph.data.frame(d1[1:2], directed = F)
> plot(g1)

> d2 <- read.table("test2.txt")
> g2 <- graph.data.frame(d2[1:2], directed = F)
> plot(g2)

・グラフの和。

> g_uni <- graph.union(g1, g2)
> plot(g_uni)

・積

> g_int <- graph.intersection(g1, g2, keep.all.vertices = FALSE)
> plot(g_int)

・差

> g_dif <- graph.difference(g1, g2)
> plot(g_dif)

・せっかくなので、和のグラフに関して、どのノードがどのグラフ由来なのか、どこが重なっていたのか色分けしてみる。

#とりあえずノードのリストを作成

> col <- vertex.attributes(g_uni)$name
> names(col)<- col
> col
A B C D E F G H I J K L M
"A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M"

#g1由来のノードを青に、g2由来のノードを赤に、被っていたノードを紫にする。

> col <- replace(col, vertex.attributes(g1)$name, "blue")
> col <- replace(col, vertex.attributes(g2)$name, "red")
> col <- replace(col, vertex.attributes(g_int)$name, "purple")
> col