【R】wikipathwayでエンリッチ解析してみたい。

・お題：WikiPathwaysという生物学的パスウェイのデータベースがあるらしい。エンリッチ解析やパスウェイの図示もできるらしい。RからCytoscapeを介して使ってみたい。

・なお、私は素人なので、ちゃんとしたことはちゃんとした人に聞いていただきたい。

・WikiPathwaysに関するRパッケージは以下のrWikiPathwaysを参照。インストールは割愛。bioconductor.org

・Cytoscapeに関しては以下を参照。ネットワーク解析に使うソフトウェア。

cytoscape.org

・今回は、気になる遺伝子のリストを得た状態からスタートし、どういったパスウェイにそれがエンリッチしているのか示唆を得て、パスウェイを図示する感じにしたい。

・まず、Cytoscapeのソフトウェアを立ち上げる。Appsから、WikiPathwaysをインストールしておく。

・次に、Rの方でライブラリを読み込む。

> library(rWikiPathways) #wikipathwayのパッケージ

> library(clusterProfiler) #エンリッチメント解析のパッケージ
> library(RCy3) #cytoscapeのパッケージ

> library(tidyverse) #Rを便利に使うパッケージ

・解析対象にできる生物種は以下のコマンドで確認できる。今回はSaccharomyces cerevisiaeでやってみる。

> get_wp_organisms()

[1] "Sus scrofa" "Solanum lycopersicum" "Saccharomyces cerevisiae" "Rattus norvegicus" "Populus trichocarpa"
[6] "Pan troglodytes" "Mus musculus" "Homo sapiens" "Gallus gallus" "Equus caballus" [11] "Drosophila melanogaster" "Danio rerio" "Canis familiaris" "Caenorhabditis elegans" "Bos taurus"
[16] "Arabidopsis thaliana"

・データセットを作る。スコアは適当に作った。

> Score <- c(1, 2, 3, 4, 5)
> Ensemble <- c("YOR317W", "YIL009W", "YLR284C", "YKR009C", "YER015W")
> TargetGene_df <-
+ data.frame(Ensemble = Ensemble,
+ Score = Score)
> TargetGene_df
Ensemble Score
1 YOR317W 1
2 YIL009W 2
3 YLR284C 3
4 YKR009C 4
5 YER015W 5

・エンリッチメント解析してみる。

> ScEnrichWp <- enrichWP(TargetGene_df$Ensemble,
+ organism = "Saccharomyces cerevisiae",
+ pAdjustMethod = "fdr",
+ pvalueCutoff = 0.1)

・なんだか怒られた。遺伝子IDの形式がパッケージで想定されているものと違っているらしい。このIDは、EntrezIDっぽいので、EntrezIDに変換してから放り込むことにする。

・ID変換には、org.Sc.sgd.dbパッケージを使う。正しいことは以下を参照。

bioconductor.org

・変換してみる。

> library("org.Sc.sgd.db") #パッケージを読み込む。

・どんな変換が可能か確認してみる。

> keytypes(org.Sc.sgd.db)
[1] "ALIAS" "COMMON" "DESCRIPTION" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS" "ENTREZID" "ENZYME" "EVIDENCE"
[10] "EVIDENCEALL" "GENENAME" "GO" "GOALL" "INTERPRO" "ONTOLOGY" "ONTOLOGYALL" "ORF" "PATH"
[19] "PFAM" "PMID" "REFSEQ" "SGD" "SMART" "UNIPROT"

・変換してみる。今回は、元の遺伝子IDの形式がENSEMBL、変換先がENTREZID。

> TargetGene_df$Entrez <-
+ mapIds(x = org.Sc.sgd.db,
+ keys = TargetGene_df$Ensemble,
+ column = "ENTREZID",
+ keytype = "ENSEMBL")
'select()' returned 1:1 mapping between keys and columns

> TargetGene_df

・それっぽいID列として、Entrezができたので、これでenrichWPしてみる。

> ScEnrichWp <- enrichWP(TargetGene_df$Entrez,
+ organism = "Saccharomyces cerevisiae",
+ pAdjustMethod = "fdr",
+ pvalueCutoff = 0.1)

・結果を見てみる。結果から逆算して遺伝子をピックアップしたので、一つの経路にガッツリenrichした。。普通はこんな結果にならない筈。。

> dotplot(ScEnrichWp, showCategory = 20)

> ScEnrichWp_df <- data.frame(ScEnrichWp)
> ScEnrichWp_df

・今回はWP91というパスウェイにエンリッチしたので、このパスウェイをcytoscape上に引っ張ってきて図示してみたい。先にcytoscapeを立ち上げている筈だけれど、一応RからCytoscapeを操作できる状態になっているのか確認するために、pingを飛ばしておく。

> cytoscapePing()
You are connected to Cytoscape!

・繋がっていることが確認できたので、cytoscapeにWP91のパスウェイをダウンロードすれば良い。一発でWP91としても良いけれど、Rからwikipathwaysのパスウェイを検索するやり方をメモしておく。

> Target_pathways <- findPathwaysByText("Fatty acid oxidation")
> Target_pathways

・パスウェイがざっと700個弱引っ掛かってきてしまった。例えば、speciesを酵母にして絞ってみる。

> Target_Scerevisiae_pathways <-
+ Target_pathways %>%
+ dplyr::filter(.,species == "Saccharomyces cerevisiae")

> Target_Scerevisiae_pathways

・35個に絞ることができた。ここから気になるパスウェイを指定しても良い。今回はエンリッチメント解析でもトップに来ていた（一つしかなかったが）WP91を選択して、パスウェイを図示する。

> Target_Scerevisiae_pathways_IDs <- Target_Scerevisiae_pathways$id
> commandsRun(paste0('wikipathways import-as-pathway id=', Target_Scerevisiae_pathways_IDs[ 1 ] ) )

・先述のコマンドをR上で実行すると、Cytoscape上でパスウェイが描かれる。

・Scoreで色を着けたいので、元のTableをCytoscapeのノード情報にくっつける。

R上のコマンド

> loadTableData(data = TargetGene_df,
+ data.key.column = "Ensemble", #最初のデータセットの時点でEnsemblの綴りを間違っていたらしい…
+ table.key.column = "Ensembl")
[1] "Success: Data loaded in defaultnode table"

Cytoscape上のノード情報