【R】パスウェイ解析を学んでみたい①

・お題：先日の遺伝子発現の解析で、発現が減った・増えた遺伝子群が出てきた。その生物学的解釈を考えるうえで、どのようなパスウェイが変わりがちか調べる解析のことを、パスウェイ解析というらしい。やってみたい。

・私はwet人間なので、パスウェイ解析も全く分からない。そこで、ネットで調べたチュートリアルをなぞって学んでみたい。今回参考にしたのは、以下のサイト。

https://bioconductor.riken.jp/packages/3.3/bioc/vignettes/ReactomePA/inst/doc/ReactomePA.html#pathway-analysis-of-ngs-data

・パスウェイは独自に構築するものではなく、すでに知られているものを利用していくことになる。一般的に用いられるパスウェイは、KEGGとREACTOMEというやつらしい。正しいことは公式サイトをご確認いただきたい。

・KEGG：

www.genome.jp

・Reactome：

reactome.org

・参考サイトでは、Reactomeを使ってパスウェイ解析をしていたので、今回はReactomeを使う。Rを介してReactomeでパスウェイ解析する場合、ReactomePA（Reactome Pathway Analysis）というパッケージを用いるらしい。

bioconductor.org

・とりあえず、ReactomePAパッケージをインストールする。このインストール、私の環境では少し時間がかかった。

> if (!require("BiocManager", quietly = TRUE))
+ install.packages("BiocManager")

・次に、解析対象にするデータセットのサンプルを読み込む。

> library(DOSE)

> data(geneList)#DOSEパッケージに入っているgeneListオブジェクトを呼び出す。

・ちなみに、このgeneListが気になる場合は以下を参照。

bioc.ism.ac.jp

・このgeneListに何が入っているかというと、Entrez Gene IDとfoldchange（たぶん）がずらっと入っているらしい。ReactomePAの解析はEntrez Gene IDが情報として必要なので、必要に応じて遺伝子のIDを変換してやる必要がある（今回はデータセットに入っているので不要）。

> library(tidyverse)

> geneList %>% head()

・次に、今回解析対象にする遺伝子セットを作成する。

> de <- names(geneList)[abs(geneList) > 1.5] #foldchangeの絶対値が1.5以上のものをdeとする。
> de %>% head()
[1] "4312" "8318" "10874" "55143" "55388" "991"

・この遺伝子IDのセットが、なんのパスウェイに集まっているか調べたい。ReactomePAパッケージのenrichPathway関数を利用する。この関数は、遺伝子のベクトルが与えられると、カットオフ（FDR）で濃縮されたパスウェイを返してくれるらしい。

> x <- enrichPathway(gene=de, pvalueCutoff=0.05, readable=T)

> x %>% summary() %>% head()

・ちっちゃくて分かりにくいが、生物学的プロセスなどが記載してあるっぽい。もうちょっと分かりやすく図示する。

> barplot(x, showCategory=10)

・調整済p値の昇順に、パスウェイが並んでいるらしい。なんだか見やすい。

・dotplotで表現することもできる。
> dotplot(x, showCategory=10)

・次に、エンリッチメントマップを描いてみる、、と、エラーが出てきた。私のような素人はエラーが出ると詰む。。要するに、enrichMapなんて関数見つからないよ！ということらしい。

> enrichMap(x, layout=igraph::layout.kamada.kawai, vertex.label.cex = 1)
Error in enrichMap(x, layout = igraph::layout.kamada.kawai, vertex.label.cex = 1) :
could not find function "enrichMap"

・調べてみると、enrichMapはReactomePAパッケージではなく、DOSEパッケージの一部だそう。DOSEパッケージはさっき読み込んだ筈なのに、、、と思ってさらに調べてみたら、私の使ったバージョンのDOSEパッケージにはenrichMap関数がないらしい。enrichMapはいったん諦める。

・パスウェイと遺伝子の関係を可視化すべく、cnetplotを描いてみる。

> cnetplot(x, categorySize="pvalue",

+ foldChange=geneList)