【R】STRINGデータを解析してみたい

・お題：タンパク質間相互作用ネットワークで、STRINGというのがある。ここからデータをとってきて、グラフを作成し、タンパク質の中心性を算出する。また、データベースから遺伝子の生存必須性データをとってきて、ID変換し、先のタンパク質-中心性の表とくっつけて、生存必須性との関係を見てみたい。

・やることは、以前やった中心性解析にデータベースからのデータ取得やID変換などを組み合わせた感じ。今回は、代表的な病原性細菌である黄色ブドウ球菌（Staphylococcus aureus）を例にやってみる。

・とりあえず、データを落とす。

> library(tidyverse)#とりあえず読み込む
> library(STRINGdb)#STRINGからデータを引っ張ってくるライブラリ

> string.db <- STRINGdb$new(version="11.5",
+ species=93061,
+ score_threshold=400,
+ network_type="full",
+ input_directory="")
> string.db
*********** STRING - https://string-db.org ***********
(Search Tool for the Retrieval of Interacting Genes/Proteins)
version: 11.5
species: 93061
............please wait............
proteins: 2892
interactions: 118610

・score_thresholdやnetwork_typeによって取れてくるデータセットは全然変わってくると思うけれど、今回は2892種類のタンパク質と118610種類の相互作用が含まれるネットワークが取れたみたい。

・ここからグラフを描くための情報を取り出す。

> full.graph <- string.db$get_graph()
> full.graph

・次にグラフを描いて少し弄る。コードは九工大の竹本先生のHPを参考にさせて頂いた。

sites.google.com

> library(igraph)#グラフをハンドリングするためにigraphパッケージを読み込む
> full.graph <- simplify(full.graph,
+ remove.multiple=T,#重複したエッジを除去
+ remove.loops=T)#自己ループを除去

・連結成分を見てみる。

> cls <- clusters(full.graph,"weak")#連結成分求める
> cls$csize
[1] 2483 2 2 2 4 2 2 2 2
> cls$no
[1] 9

・9つの島があるらしい。そのうち一つが巨大であり、2483個のタンパク質が含まれている。解析対象はこの巨大連結成分だけにする。

> sub.graph <- induced_subgraph(full.graph,
+ V(full.graph)$name[cls$membership == which(cls$csize == max(cls$csize))[ [1] ] ])

・次に、中心性をいろいろを求めてみる。

> degree <- degree(g)
> eigen <- evcent(g)$vector
> page <- page.rank(g)$vector
> closen <- closeness(g)
> between <- betweenness(g)
> subgraph <- subgraph.centrality(g)

> nprop <- data.frame(V(g)$name, degree, eigen, page, closen, between, subgraph)
> nprop <- nprop %>%
+ dplyr::rename("STRING_id" = V.g..name)

> nprop %>% head(.)

> nprop %>% dim(.)
[1] 2483 7
・これで、データベースからデータをダウンロードして、PPIネットワークを作成し、中心性を算出するところまでいった。

・次に、生存必須性データをダウンロードする。今回は以下のデータベースからデータを得た。

v3.ogee.info

・csvをダウンロードして、カレントディレクトリに保存し、Rに読み込む。

> Sa_OGEE <-
+ read.csv("Saureus_essgenes.csv",
+ stringsAsFactors = FALSE)
> Sa_OGEE %>% head(.)

・この中で使う情報は、locusとessentialityのみ。essentialityの様子を見てみる。

> Sa_OGEE$essentiality %>% table(.)
DE E NE
190 671 4751

・DEはEssentialでもNon-Essentialでもないっぽいので、除外する。

> Sa_OGEE <- Sa_OGEE %>%
+ dplyr::select(., c(locus, essentiality)) %>%
+ dplyr::filter(., essentiality == "E"|essentiality == "NE")
> Sa_OGEE %>% head(.,10)

> Sa_OGEE %>% dim()
[1] 5422 2

・STRINGのPPIネットワークのノード数よりも、生存必須性データの方が多い。。

・生存必須性の情報とSTRINGのID情報が違うので、IDを変換する。STRINGからaliasデータを落として、変換リストを作成する。

> initial_alias_df <- read.table(paste0(tempdir(), "/", "93061.protein.aliases.v11.5.txt.gz"),
+ skip = 1, sep = "\t", header = FALSE, quote = "",
+ stringsAsFactors = FALSE, fill = TRUE,
+ col.names = c("STRING_id", "locus", "sources"))

> string_locus <- initial_alias_df %>% dplyr::filter(., sources == "RefSeq_locus")
> string_locus

・変換リストができたので、くっつけていく。中心性のデータ（nprop、IDはSTRING_id）に、変換リスト（string_locus、IDはSTRING_idとlocus）をくっつけ、さらに生存必須性のデータ（Sa_OGEE、IDはlocus）をくっつける。今回は中心性データが少ないので、中心性データに合わせる形にする。

> nprop <- nprop %>%
+ dplyr::left_join(., string_locus, "STRING_id") %>%
+ dplyr::left_join(., Sa_OGEE, "locus")

・このテーブルを使って、生存必須性の有無と中心性の関係を、中心性の種類ごとにBoxplotで示してみる。

> library(patchwork)
> g <- ggplot(nprop)
> 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)

・個人的には結構傾向が出ているように感じた。せっかくなので、それぞれの中心性で閾値を設けて生存必須性の有無を判別する場合のROCカーブを描いてみる。

> library(pROC)
> roc1 <- roc(ess_score~nprop$degree, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> roc2 <- roc(ess_score~nprop$eigen, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> roc3 <- roc(ess_score~nprop$page, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> roc4 <- roc(ess_score~nprop$closen, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> roc5 <- roc(ess_score~nprop$between, ci = TRUE)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
> roc6 <- roc(ess_score~nprop$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")