・お題：先日、参考サイトをなぞって大腸菌のネットワークの中心性指標と生存必須性の関係を見てみた。その際は参考サイトで提供されているデータセットを使って実行してみたが、自分でデータセットをとって来るところからやってみたい。
・今回の参考元サイト。解析手法はこちらをなぞる。
・今回は、黄色ブドウ球菌（Staphylococcus aureus）のデータを公共データベースから落としてきて解析してみる。
・解析には、タンパク質間相互作用（PPI）ネットワーク情報と、正解ラベルになる生存必須性の情報が必要になる。PPIネットワーク情報はSTRING、生存必須性情報はOGEEから落としてくる。
・ダウンロードと解凍は割愛。Staphylococcus aureusで検索してそれっぽいデータを落としてきて、R.utilsのgunzipなどすれば解凍できる。
・ライブラリとデータを読み込む。
> library(tidyverse)
> #PPIネットワークデータ
> d <- read.table("93061.protein.links.v11.5.txt", header = TRUE)
> d %>% head()
> d %>% nrow()
[1] 621232
> #生存必須性データ
> ess<-read.csv("Staphylococcus aureus subsp. aureus NCTC 8325_genes.csv", header = TRUE)
> ess %>% head()
> ess %>% nrow()
[1] 2891
> table(ess$essentiality)
E NE
350 2541
・ネットワーク情報のタンパク質名と生存必須性情報のlocusがちょっと違う。前者はtaxaIDがくっついているので、"."でカチ割って表記を合わせる。"."自体が正規表現なので、抜ける必要がある。
> d <- separate(d, protein1, c("taxaID1", "protein1"), sep="\\.")
> d <- separate(d, protein2, c("taxaID2", "protein2"), sep="\\.")
> d %>% head()
・PPIネットワークからグラフを読み込む。
> library(igraph)
> g <- graph.data.frame(d[c(2,4)], directed=T)
> E(g)$weight <- d[ [ 5 ] ]
> g <- igraph::simplify(g,
+ remove.multiple=T,
+ remove.loops=T)
> V(g)
・最大連結成分（最も大きいネットワークの島）を探す。
> cls <- clusters(g, "weak")
> cls$csize
[1] 2598 2
> g <- induced_subgraph(g,
+ V(g)$name[cls$membership == which( cls$csize == max( cls$csize ))[ [ 1 ] ] ] )
> V(g)
・中心性指標を求める。
> # 次数中心性（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)
>
> nprop <- data.frame(V(g)$name, degree, eigen, page, closen, between, subgraph)
> names(nprop)[ [ 1 ] ] <- "locus"
> nprop %>% head()
・中心性指標のテーブルと生存必須性のテーブルをlocusを指標としてくっつける。
> d <- merge(ess, nprop, by = "locus")
> d[,7] %>% table()
E NE
345 2253
・箱ひげ図に示してみる。参考サイトと少し違う手法を使って、ggplotとpatchworkの練習も兼ねてみた。
> library(patchwork)
> g <- ggplot(d)
> g1 <- g + geom_boxplot(aes(x = essentiality, y = log(degree)))
> g2 <- g + geom_boxplot(aes(x = essentiality, y = log(eigen)))
> g3 <- g + geom_boxplot(aes(x = essentiality, y = log(page)))
> g4 <- g + geom_boxplot(aes(x = essentiality, y = log(closen)))
> g5 <- g + geom_boxplot(aes(x = essentiality, y = log(between + 1)))
> g6 <- g + geom_boxplot(aes(x = essentiality, y = log(subgraph)))
> (g1/g2)|(g3/g4)|(g5/g6)
・いずれの中心性も必須遺伝子の方が高い傾向にあるものの、近接中心性と媒介中心性はなんだかかなり被っているように見える。wilcoxonの順位和検定をしてみる。
> #必須と非必須の中心性をそれぞれ分けて表にする、
> ess_nprop <- d[d$essentiality=="E",c(10:15)]
> noness_nprop <- d[d$essentiality=="NE",c(10:15)]
> #統計
> d2 <- data.frame(matrix(ncol = ncol(ess_nprop), nrow = 0))
> colnames(d2) <- names(ess_nprop)
> for (n in 1:ncol(d2)){
+ d2[1, n] <- wilcox.test(ess_nprop[,n], noness_nprop[,n])$p.value
+ }
> #ついでに、それぞれの中央値も算出してまとめる。
> 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
・個別の中心性を指標に必須性を予測した場合のROCカーブを描いてみる。
> # 必須なら1，非必須なら0とする。
> ess_score <- ifelse(d$essentiality == "E",1,0)
> ess_score %>% head()
[1] 1 1 1 0 1 1
> # ROCカーブを書くためのパッケージを読み込む
> library(pROC)
>
> roc1 <- roc(ess_score~d$degree, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> roc2 <- roc(ess_score~d$eigen, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> roc3 <- roc(ess_score~d$page, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> roc4 <- roc(ess_score~d$closen, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> roc5 <- roc(ess_score~d$between, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> roc6 <- roc(ess_score~d$subgraph, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> layout(matrix(1:6, ncol=3))
> plot(roc1, print.auc = TRUE, legacy.axes = TRUE, main = "Degree")
> plot(roc2, print.auc = TRUE, legacy.axes = TRUE, main = "Eigenvector")
> plot(roc3, print.auc = TRUE, legacy.axes = TRUE, main = "PageRank")
> plot(roc4, print.auc = TRUE, legacy.axes = TRUE, main = "Closeness")
> plot(roc5, print.auc = TRUE, legacy.axes = TRUE, main = "Betweenness")
> plot(roc6, print.auc = TRUE, legacy.axes = TRUE, main = "Subgraph")
・近接中心性と媒介中心性のゆがみが激しい。AUCが見づらいので、個別に算出してみる。
> roc1$auc#次数中心性
Area under the curve: 0.831
> roc2$auc#固有ベクトル中心性
Area under the curve: 0.8638
> roc3$auc#PageRank
Area under the curve: 0.8286
> roc4$auc#近接中心性
Area under the curve: 0.7317
> roc5$auc#媒介中心性
Area under the curve: 0.6211
> roc6$auc#サブグラフ中心性
Area under the curve: 0.8574
・箱ひげ図のイメージ同様、媒介中心性と近接中心性が指標としてかなり怪しい結果になった。対して、固有ベクトル中心性とサブグラフ中心性がスコア的にはよかった。
・個別の遺伝子の必須性とネットワークでの位置づけに関してみていくと、発見がありそう。それぞれの中心性に対する理解が解釈には必要そう。中心性を組み合わせたり、遺伝子の特徴量として機械学習に与えることで必須性予測の性能が向上するような気もするけれど、どうなんだろうか。。
おわり。
