- clusterProfilerを用いた KEGGパスウェイ解析
- ReactomePA を用いたReactomeパスウェイ解析
この記事でご紹介するのは上記の二つです！clusterProfilerとReactomePA はどちらもRのパッケージです．
RNA-seq解析やマイクロアレイ解析では発現が増加したり，減少したりした遺伝子群（発現変動遺伝子）を抽出し，それらがどのような機能を持っているのかを調べます．
パスウェイ解析はその機能解析のうちの一つの手法であり，興味のある遺伝子群（今回は発現変動遺伝子群）がどんなパスウェイに所属しているのかを調べる解析方法です．
基本的な考え方はGene Ontology解析（GO解析）と同様です．GO解析についてはこちらが詳しいです．https://bi.biopapyrus.jp/pathway/go/
KEGG パスウェイ解析
パスウェイ解析をするならまずはこれ，というぐらい有名なデータベースがKEGG（Kyoto Encyclopedia of Genes and Genomes）です！KEGGはオープンソース，オープンアクセスのパスウェイデータベースです．
実際にパスウェイ解析を行えるようにサンプルデータを用意しましたので，よろしければ以下からダウンロードしてご使用ください．
https://github.com/koji11235/pathway_enrichment_analysis_r/blob/main/sample_gene_data.csv
サンプルデータはairwayというデータセットから調整P-valueが0.01以下の遺伝子だけを抽出したものです．
遺伝子の識別子であるEntrez Gene ID，発現変動量であるlog2FoldChangeなどが含まれています．
この後の解析はこのサンプルデータを用いて行なっています．
サンプルデータファイルをワーキングディレクトリに置いてある状況を想定しています．
#ライブラリの読み込み
library(clusterProfiler)
#サンプルデータ読み込み
sample_data<-read.csv("sample_gene_data.csv")
#サンプルデータの中のEntrez_Gene_IDを抽出
gene_IDs<-sample_data$Entrez_Gene_ID
#KEGGエンリッチメント解析を実行
kegg_enrichment_result <- enrichKEGG(gene=gene_IDs,pvalueCutoff=0.05)
#結果の表示
head(kegg_enrichment_result, n=10)
ほんの数行でパスウェイ解析を行うことができました！
結果のDescriptionの列に表示されるのがパスウェイ解析で得られたパスウェイの説明，IDの列に表示される「”hsa”＋5桁の数字」がそのパスウェイのIDです．
次に得られたパスウェイがどんなものかをパスウェイIDを使って見てみます．
p valueの最も低いFocal adhesionのパスウェイ（hsa04510）をKEGGに見に行って見ましょう！
browseKEGG(kegg_enrichment_result,"hsa04510")
上記のコードを実行するとKEGGの対象のパスウェイのページに飛ぶことができます．（結果のページ）さらにそのパスウェイの図の中で今回解析に用いた遺伝子が色付けて表示してくれています！
このパスウェイの図をダウンロードしてみましょう！
ついでに遺伝子の発現変動の量も描画しちゃいましょう！
pathviewというパッケージを使用します．
#ライブラリの読み込み
library(pathview)
#サンプルデータからlog2FoldChangeを抽出
logFC <- sample_data$log2FoldChange
#logFCとEntrez_Gene_IDを対応させる
names(logFC) <- sample_data$Entrez_Gene_ID
#パスウェイの画像のダウンロード
pathview(gene.data = logFC, pathway.id = "hsa04510", limit = list(gene=2, cpd=1))
上記のコードを実行すると下のような画像が現在のディレクトリにダウンロードされます！
ダウンロードされるファイル名は”hsa04510.pathview.png”となっています．
色はlog2FoldChangeを表しています．
パスウェイ中でどの遺伝子の発現が変化しているのかが一目でわかりますね！
KEGGパスウェイ解析については以上になります．
Reactome パスウェイ解析
次はReactomeパスウェイ解析です！ReactomeはKEGGと同様にオープンソース，オープンアクセスのパスウェイデータベースです．モダンなデザインのデータベースで，使用しているイラストも非常に綺麗なのが特徴です．
KEGGパスウェイ解析で行なったのと同様にサンプルデータを用いて解析を行います．
ReactomePAというパッケージを使用します．
#ライブラリの読み込み
library(ReactomePA)
#サンプルデータの読み込み
sample_data<-read.csv("sample_gene_data.csv")
#サンプルデータの中のEntrez_Gene_IDを抽出
gene_IDs<-sample_data$Entrez_Gene_ID
#サンプルデータからlog2FoldChangeを抽出
logFC <- sample_data$log2FoldChange
#logFCとEntrez_Gene_IDを対応させる
names(logFC) <- sample_data$Entrez_Gene_ID
#Reactomeパスウェイ解析を実行
Reactome_enrichment_result <- enrichPathway(gene=gene_IDs,pvalueCutoff=0.05, readable=T)
#結果の表示
head(as.data.frame(Reactome_enrichment_result))
KEGGパスウェイ解析と違うのはenrichPathway()の部分のみです．
こちらもほんの数行でパスウェイ解析を実行できました！
結果のDescriptionの列に表示されるのがパスウェイ解析で得られたパスウェイの説明，IDの列に表示される「”R-HSA”＋7桁の数字」がそのパスウェイのIDです．
ReactomePAがすごいのはここからで，様々な種類の可視化に対応しています．
#バープロット
barplot(Reactome_enrichment_result, showCategory=8,x = "Count")
最もよく使われるパスウェイ解析の結果の可視化がこのバープロットです．棒グラフの高さは遺伝子の数を表しています．x＝’GeneRatio’ とすると「そのパスウェイに含まれる全遺伝子の中で今回の解析で発現変動していた遺伝子の数」を割合として描画してくれます．
#ドットプロット
dotplot(Reactome_enrichment_result, showCategory=15)
ドットプロットは横軸でパスウェイに含まれている全遺伝子の中で今回の解析で発現変動していた遺伝子の割合を，
ドットの大きさでパスウェイ中で今回の解析で発現変動していた遺伝子の数を表しています．
先ほどのバープロットよりも表現できる情報が一つ増えるのがこの可視化を用いるメリットです．
#エンリッチメントマッププロット
emapplot(Reactome_enrichment_result)
エンリッチメントマップは各パスウェイを頂点にしたネットワーク図です．パスウェイ同士の繋がりはパスウェイに含まれていた遺伝子の重複度合いを基準にして書かれています．つまり，繋がっているパスウェイ同士は多くの共通の遺伝子を持つということです．
同じ遺伝子を共有する結果として似た機能を持ったパスウェイ同士がクラスターを形成するためパスウェイ同士の関係性が非常に見やすいです．
#cnetプロット
cnetplot(Reactome_enrichment_result, categorySize="pvalue", foldChange=logFC,showCategory = 8)
cnetプロットはパスウェイの名前を中心にしてその周りにパスウェイに所属する遺伝子が繋がっているネットワーク図です．
パスウェイが同じ遺伝子を共有していた場合は一つの遺伝子から複数のパスウェイに繋がります．上のエンリッチメントマップのパスウェイ間の繋がりの中身まで描画したプロットと言えます．こちらの方が表現できる情報量はかなり多いです．
日本語訳が見つからなかったのですが，Concept Network Plotの略かと思います
．
ちなみに上記の4つの可視化ですが，clusterProfilerで行なったKEGGの結果もこれらの関数で表示することができます．これはclusterProfilerとReactomePAの作者が同じだからですね．作者のGuangchuang Yuさんに感謝です．
もう一つちなみに上記4つの可視化は全てggplotパッケージによって描画されていますのでggplotを用いたより細かな調整を行うことができます．またemapplotとcnetplotはネットワーク描画に特化したggraph パッケージによって描画されています．これら2つについてはggraph を用いたノードやエッジの追加を行うことができます．
ReactomeでKEGGのようなパスウェイの図をダウンロードする
ReactomeのAPIでDiagram exporterを用いれば可能です！
Reactomeパスウェイ解析の結果得られた”Signaling by Receptor Tyrosine Kinases”（R-HSA-9006934）の図を取得してみましょう．
ブラウザに以下のURLを入力するとR-HSA-900693のパスウェイの図がダウンロードできます．
https://reactome.org/ContentService/exporter/diagram/R-HSA-9006934.png
めちゃくちゃ綺麗な図が得られましたね！
“https://reactome.org/ContentService/exporter/diagram/”+”パスウェイID”+”拡張子”
とするとパスウェイIDに対応した図を得ることができます．
このダウンロードをRでやってみます．
#パスウェイIDを取得（R-HSA-9006934）
pathwayID<-Reactome_enrichment_result[1]$ID
#URLを設定
URL<-paste0("https://reactome.org/ContentService/exporter/diagram/",pathwayID,".png")
#保存するファイル名を設定
output_filename<-paste0(pathwayID,".png")
#ダウンロードを実行
download.file(URL,output_filename)
上のコードを実行すると作業ディレクトリに”R-HSA-9006934.png”というファイル名でダウンロードできます．
KEGGでやったのと同じようにこれにどの遺伝子が含まれていたのかを表示させてみましょう．URLに ?flg=GeneName,GeneName という形で遺伝子名を渡せば該当する遺伝子をハイライトした図を表示してくれます．
#文字列加工用のライブラリー
library(stringr)
#パスウェイIDを取得（R-HSA-9006934）
pathwayID<-Reactome_enrichment_result[1]$ID
#含まれている遺伝子を抽出，加工
genelist<-str_split(Reactome_enrichment_result[1]$geneID,"/")[[1]]
genelist.str<-str_c(genelist,collapse = ",")
#URLを設定
URL<-paste0("https://reactome.org/ContentService/exporter/diagram/",pathwayID,".png","?flg=",genelist.str)
#保存するファイル名を設定
output_filename<-paste0(pathwayID,"_decorated",".png")
#ダウンロードを実行
download.file(URL,output_filename)
上のような画像がダウンロードできます！
上記のコードで作成したURLが以下です．
おわりに
ご覧いただきましてありがとうございました！少しでもお役に立つ情報であれば幸いです．他の手法でパスウェイ解析を行う方法がありましたらぜひご教示ください．
参考ページ
https://www.genome.jp/kegg/kegg_ja.html
https://bi.biopapyrus.jp/pathway/go/
https://bioconductor.org/packages/release/bioc/vignettes/ReactomePA/inst/doc/ReactomePA.html
https://bioinformatics-core-shared-training.github.io/cruk-summer-school-2018/RNASeq2018/html/06_Gene_set_testing.nb.html
https://reactome.org/dev/content-service/diagram-exporter
https://bioc.ism.ac.jp/packages/3.7/bioc/vignettes/enrichplot/inst/doc/enrichplot.html#enrichment-map
Citation
Guangchuang Yu, Li-Gen Wang, Yanyan Han and Qing-Yu He.
clusterProfiler: an R package for comparing biological
themes among gene clusters. OMICS: A Journal of
Integrative Biology 2012, 16(5):284-287
Guangchuang Yu, Qing-Yu He. ReactomePA: an
R/Bioconductor package for reactome pathway analysis and
visualization. Molecular BioSystems 2016, 12(2):477-479