受託コンサルティングチームの kutsuy でございます。 前編では ORA（Over-Representation Analysis）と GSEA（Gene Set Enrichment Analysis）の事前準備をしました。 さっそく ORA の例として、enrichGO による解析を実行してみましょう。
動作環境（前編分を含む）
- R version 4.3.0
- airway v1.20.0
- DESeq2 v1.40.2
- clusterProfiler v4.8.3
- org.Hs.eg.db v3.17.0
- dplyr v1.1.4
- enrichplot v1.20.3
ORA（Over-Representation Analysis）
# 有意な遺伝子を抽出（例：padj < 0.05） > sig_genes <- res2 %>% filter(padj < 0.05) %>% pull(ENTREZID) > length(sig_genes) [1] 2541 > head(sig_genes) [1] "3927" "90293" "8837" "9108" "381" "23028" > ego <- enrichGO(gene = sig_genes, + OrgDb = org.Hs.eg.db, + keyType = "ENTREZID", + ont = "BP", + pAdjustMethod = "BH", + pvalueCutoff = 0.05, + qvalueCutoff = 0.2, + readable = TRUE)
発現変動が有意（padj < 0.05）な 2541 個の Entrez ID のベクトルを入力としています。 試しに、ego に格納された結果の一部を表示してみます（読み飛ばしていただいても構いません）。
> class(ego) [1] "enrichResult" attr(,"package") [1] "DOSE" > ego # # over-representation test # #...@organism Homo sapiens #...@ontology BP #...@keytype ENTREZID #...@gene chr [1:2540] "3927" "90293" "8837" "9108" "381" "23028" "5965" "5166" "340273" "7982" "93655" "5976" "56919" ... #...pvalues adjusted by 'BH' with cutoff <0.05 #...692 enriched terms found 'data.frame': 692 obs. of 9 variables: $ ID : chr "GO:0030198" "GO:0043062" "GO:0045229" "GO:0061448" ... $ Description: chr "extracellular matrix organization" "extracellular structure organization" "external encapsulating structure organization" "connective tissue development" ... $ GeneRatio : chr "83/2275" "83/2275" "83/2275" "75/2275" ... $ BgRatio : chr "314/18614" "315/18614" "317/18614" "275/18614" ... $ pvalue : num 4.05e-12 4.85e-12 6.94e-12 8.57e-12 9.90e-11 ... $ p.adjust : num 1.32e-08 1.32e-08 1.32e-08 1.32e-08 1.22e-07 ... $ qvalue : num 1.04e-08 1.04e-08 1.04e-08 1.04e-08 9.59e-08 ... $ geneID : chr "CFLAR/ST7/PHLDB1/TNFRSF1B/GPM6B/COL9A2/COL11A1/BCL3/ITGA8/COL5A3/COL4A4/B4GALT1/ADAMTS2/LAMB1/TGFB2/CST3/RGCC/M"| __truncated__ "CFLAR/ST7/PHLDB1/TNFRSF1B/GPM6B/COL9A2/COL11A1/BCL3/ITGA8/COL5A3/COL4A4/B4GALT1/ADAMTS2/LAMB1/TGFB2/CST3/RGCC/M"| __truncated__ "CFLAR/ST7/PHLDB1/TNFRSF1B/GPM6B/COL9A2/COL11A1/BCL3/ITGA8/COL5A3/COL4A4/B4GALT1/ADAMTS2/LAMB1/TGFB2/CST3/RGCC/M"| __truncated__ "CFLAR/PKD1/CHRDL2/COL11A1/RASAL2/OXCT1/HIF1A/MAPK3/BPNT2/CCN4/DYRK1B/COMP/ITGB8/HOXA3/GLI3/ACTA2/COL1A1/ZBTB16/"| __truncated__ ... $ Count : int 83 83 83 75 102 99 106 79 93 84 ... #...Citation T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
次に、GSEA の例として、gseGO による解析を実行してみましょう。
GSEA（Gene Set Enrichment Analysis）
# fold change に基づくランキング
> geneList <- res2$log2FoldChange
> names(geneList) <- res2$ENTREZID
> length(geneList)
[1] 37362
> geneList <- sort(geneList, decreasing = TRUE)
> length(geneList)
[1] 22661
> head(geneList)
247 7704 3000 10218 401613 79689
9.609736 7.176082 6.026012 5.796834 5.535285 5.287032
> gse <- gseGO(geneList = geneList,
+ OrgDb = org.Hs.eg.db,
+ keyType = "ENTREZID",
+ ont = "BP",
+ minGSSize = 10,
+ maxGSSize = 500,
+ pvalueCutoff = 0.05,
+ verbose = FALSE)
There were 12 warnings (use warnings() to see them)
logFC の大きい順に sort 関数で並び替えたときに、値が NA の遺伝子IDが除外されて 22661 個となりました。 logFC が大きい順に並んだ 22661 個の logFC値（Entrez ID を名前として持つ）のベクトルを入力としています。 試しに、gse に格納された結果の一部を表示してみます（読み飛ばしていただいても構いません）。
> class(gse) [1] "gseaResult" attr(,"package") [1] "DOSE" > gse # # Gene Set Enrichment Analysis # #...@organism Homo sapiens #...@setType BP #...@keytype ENTREZID #...@geneList Named num [1:22661] 9.61 7.18 6.03 5.8 5.54 ... - attr(*, "names")= chr [1:22661] "247" "7704" "3000" "10218" ... #...nPerm #...pvalues adjusted by 'BH' with cutoff <0.05 #...11 enriched terms found 'data.frame': 11 obs. of 11 variables: $ ID : chr "GO:0006346" "GO:0140718" "GO:0071294" "GO:0035357" ... $ Description : chr "DNA methylation-dependent heterochromatin formation" "facultative heterochromatin formation" "cellular response to zinc ion" "peroxisome proliferator activated receptor signaling pathway" ... $ setSize : int 30 43 18 23 23 14 37 46 111 236 ... $ enrichmentScore: num 0.778 0.699 0.815 0.773 0.767 ... $ NES : num 2.26 2.18 2.07 2.12 2.1 ... $ pvalue : num 1.04e-07 1.22e-06 6.82e-06 1.42e-05 1.97e-05 ... $ p.adjust : num 0.000638 0.003739 0.013894 0.021712 0.024115 ... $ qvalue : num 0.000626 0.003669 0.013634 0.021306 0.023664 ... $ rank : int 131 131 1074 1605 1074 1002 1074 1467 3679 4002 ... $ leading_edge : chr "tags=17%, list=1%, signal=17%" "tags=12%, list=1%, signal=12%" "tags=44%, list=5%, signal=42%" "tags=30%, list=7%, signal=28%" ... $ core_enrichment: chr "125997/729458/284428/653657/653656" "125997/729458/284428/653657/653656" "4501/4496/4495/4499/4502/4493/3777/4494" "247/3952/440503/59344/10370/55885/5468" ... #...Citation T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
最後にそれぞれの解析結果を可視化してみましょう。 enrichplot パッケージを使用しますが、こちらについては別記事でもご紹介しております。 enrichGO の例としてドットプロットを作図してみます。
ORAの可視化
> library(enrichplot) # ドットプロット > dotplot(ego, showCategory = 20)
他にも1行のコマンドでいろんな作図ができますのでぜひお試しください。 次に、gseGO の可視化として、GSEA特有の2つのプロットを作図してみます。
GSEAの可視化
# ridgeplot（GSEA特有のプロット） > ridgeplot(gse) Picking joint bandwidth of 0.381
上図はリッジプロットです。 有意と判定された11個のGOタームが縦軸に並んでいます（デフォルトは上位30個）。 横軸は遺伝子の log2 fold change（またはスコア）です。 また、各 GO term に関連する遺伝子のスコア分布を密度曲線として描画します。 例えば曲線が左に偏っているタームは、関連遺伝子が下降している（log2FC < 0）となるので、 「trt（処理有り）で発現が低下」している遺伝子を多く含む機能と解釈されます。 プロットの色も補正 p 値を反映していてわかりやすいですね。 見方が特殊で混乱しそうですが、落ち着いて解釈するようにしましょう。
また、別の関数でGSEA プロットを作図してみます。 こちらは興味のある個別のGOタームについての作図になります。 最も統計値が低く有意にエンリッチしていたタームを指定して実行してみましょう。
# GSEAプロファイル（個別のGO termを指定） gseaplot2(gse, geneSetID = gse$ID[1], title = gse$Description[1])
こちらもGSEA独特の図になっています。 上部の黄緑色の折れ線グラフは Running Enrichment Score を表します。 頂点が左寄りなら上昇遺伝子に多く含まれる、右寄りなら下降遺伝子に多く含まれることを意味します。 中央の黒色の縦線は遺伝子リストの中で、対象の GO term に含まれる遺伝子の位置を示します。 下部のカラーグラデーションは遺伝子リストの fold change（またはスコア）を色で表現しています。 通常、赤系が上昇、青系が下降を示します。
出力した図では明らかにピークが左側に固まっているため、 指定したGOタームでは発現上昇した遺伝子が多く（？）含まれることがわかります。
ここでリッジプロットとGSEAプロットの解釈に矛盾が生じていると思われた方は鋭いです。
これはそれぞれのプロットの視点が異なることが原因だと考えられます。
gseaplot2() は特定の GO term に注目して、その term に含まれる遺伝子が fold change 順にどこに分布しているかを示します。
ridgeplot() は複数の GO term を並列に比較し、それぞれの term に含まれる遺伝子の fold change 分布を密度曲線で示します。但し、密度曲線の形状は fold change の分布に依存するため、少数の強い上昇遺伝子がいても、全体が下降傾向なら左に偏ることがあります。
おそらく状況としては、「GO term に含まれる遺伝子のうち、一部が強く上昇している。しかし、全体としては fold change が負に偏っている」ということが起こっている可能性があります。 一般的には gseaplot2() を優先して解釈するようですが、そのようなデータであるということは解析者が認識しておくべきだと考えます。
とっても話が長くなりましたが、お役に立てそうな情報はありましたでしょうか。 GSEA解析については分かりやすく説明している解説があまり見当たらないので、今回調査も含めてまとめてみた次第です。 皆様のお困りごとが少しでも解決されれば幸いでございます。
どうでもいいですが、よくGSEAをGESAと間違えます。 それでは皆様ごきげんよう。
