・お題：ネットワークの大事なところを探すことを、中心性解析と呼ぶらしい。生物の分子機能のネットワークの中心性解析を使うと、生体維持に大切なところを見つけられたりするみたい。ぜひやってみたい。
・先日ネットワーク解析を学びたいというシリーズの記事を書いた際に、参考サイトの著者でもいらっしゃる九工大の竹本先生から教えて頂いた資料を参考に学んでみる。大変感謝です。
・今回の参考元サイト。正しいことはぜひそちらをご確認いただきたい。
・生物ネットワークに関して学びたい方は、竹本先生の貴著がおすすめ。
・まず今回使用するデータセットをダウンロードする。竹本先生のHP（以下）にて、今回使用するデータセットをダウンロードし、zipを解凍し、カレントディレクトリにecoli_ppi_Hu_etal_2009.txtとecoli_proteins_essentiality_Baba2006MSB.txtを保存した。
・ライブラリを立ち上げる。
> library(igraph)
・次に、データを読み込む。ネットワークとして読み込むデータはtextファイルで、ノードとノードの繋がり（エッジ）がずらっと入ったリストになっている。
> d <- read.table("ecoli_ppi_Hu_etal_2009.txt")
> g <- simplify(graph.data.frame(d, directed=F),#無向グラフとして読み込む
+ remove.multiple=T,#重複したエッジを除去
+ remove.loops=T)#自己ループを除去
> V(g)
・1757ノードあるので、図示するのは厳しそう。。最大連結成分（最も大きな島）を解析対象とする。まず連結成分を求める。
> cls <- clusters(g,"weak")
>cls
中略
・37の連結成分があり、さらに1673ノードが含まれている圧倒的に巨大な連結成分がある。この最大連結成分を解析対象にする。最大サイズを求め、最大サイズの連結成分を取得し、サブグラフとして抽出する。
> g <- induced_subgraph(g, V(g)$name[ cls$membership==which(cls$csize==max( cls$csize))[ [1] ] ] )
> V(g)
・gのノード数が1673になった。
・中心性を求める。
> # 次数中心性（degree centrality）
> degree <- degree(g)
> # 固有ベクトル中心性（eigenvector centraliry）
> eigen <- evcent(g)$vector
> # PageRank
> page <- page.rank(g)$vector
> # 近接中心性（closeness centrality）
> closen <- closeness(g)
> # 媒介中心性（betweenness centrality）
> between <- betweenness(g)
> # サブグラフ中心性（subgraph centrality）
> subgraph <- subgraph.centrality(g)
・一つ取り出して中身を見てみる。
> page
・ひとつひとつのノード（遺伝子）の中心性がずらっと並んでいる。この後、各中心性と生存必須性の関係を調べていくことになるのだけれど、個別に調べるのも大変なので、データフレームにまとめる。
> nprop <- data.frame(V(g)$name, degree, eigen, page, closen, between, subgraph)
> names(nprop)[ [1] ] <- "gene"#最初の列の名前をgeneに変えておく。生存必須性のtableの項目名と合わせる。
> nprop
・必須性のデータを読み込む。
> ess <- read.table("ecoli_proteins_essentiality_Baba2006MSB.txt", header = T)
> ess
gene essential
1 b0001 N
2 b0002 N
3 b0003 N
4 b0004 N
…
・必須（E），非必須（N），不明（u）とのこと。内訳を覗いてみる。
> colSums(table(ess))
E N u
234 3108 1019
> round(colSums(table(ess))/sum(table(ess)),2)
E N u
0.05 0.71 0.23
・思ったよりもuが多い印象。。次に、中心性指標のデータと生存必須性のデータをマージする。
> d <- merge(ess, nprop,by="gene")
> d
・この後essentialと中心性指標の関係を調べていくことになるけれど、essentialが不明のところは判断できないので、削除する。
> d <- d[d$essential!="u",]
> d
・縦軸にlog(中心性)、横軸に生存必須性（E/N）をとって、Boxplotを描いてみる。今回ggplotの練習も兼ねて実行したいので、tidyverseとpatchworkをアクチベートしておく。
> library(tidyverse)
> library(patchwork)
・対数をとる前に、中心性の最小値を確認しておく（ゼロ以下の値が含まれないか見ておく）。
> apply(d, 2, min)#2番目の引数が2だと列。
・geneとessentialはさておき、betweenは最小値がゼロになっており、対数をとれない。１足してから対数をとることにして、グラフを描く。
> g <- ggplot(d)
> g1 <- g + geom_boxplot(aes(x = essential, y = log(degree)))
> g2 <- g + geom_boxplot(aes(x = essential, y = log(eigen)))
> g3 <- g + geom_boxplot(aes(x = essential, y = log(page)))
> g4 <- g + geom_boxplot(aes(x = essential, y = log(closen)))
> g5 <- g + geom_boxplot(aes(x = essential, y = log(between + 1)))
> g6 <- g + geom_boxplot(aes(x = essential, y = log(subgraph)))
> g1+g2+g3+g4+g5+g6
・イメージ通り、生存必須な遺伝子ほどネットワークにおける中心性が高い雰囲気がある。統計的にどうなのか調べてみる。今回は、wilcoxonの順位和検定を用いる。以下のサイトが分かりやすかったので、気になる方はご覧ください。
・前準備として、必須タンパク質（d$essential=="E"）の中心性指標だけのdata.frameと、非必須性タンパク質（d$essential=="N"）の中心性指標だけのdata.frameを作っておく。ついでに、計算に不要な１列目と２列目を除いておく。
> ess_nprop <- d[d$essential=="E", c(-1,-2)]
> noness_nprop <- d[d$essential=="N", c(-1,-2)]
・ざっくりp値を計算してまとめる。
> d2 <- data.frame(matrix(ncol = ncol(d)-2, nrow = 0))#空データフレームを作る
> colnames(d2) <- names(d)[c(-1,-2)]#列名を入れておく
> for (n in 1:ncol(d2)){
+ d2[1, n] <- wilcox.test(ess_nprop[,n], noness_nprop[,n])$p.value
+ }
> d2
・元サイトとp値の結果が違っていて焦ったが、例えばdegreeに関してwilcoxon.testを実行した際の結果はp-value < 2.2e-16であり、ここからp.valueを抜き出すと2.440691e-30になっていた。表示している数値と内部に格納されている数値の違いっぽい。。
> wilcox.test(ess_nprop$degree, noness_nprop$degree)
Wilcoxon rank sum test with continuity correction
data: ess_nprop$degree and noness_nprop$degree
W = 171775, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0
> wilcox.test(ess_nprop$degree, noness_nprop$degree)$p.value
[1] 2.440691e-30
・せっかくなので、中央値も算出してdata.frameにまとめてみる。
> d2 <- rbind(d2,
+ apply(ess_nprop, 2, median),
+ apply(noness_nprop, 2, median))
>
> d2 <- d2 %>% t() %>% as.data.frame()
> colnames(d2) <- c("wilcox_p","E_median","N_median")
> d2
つづく。
