clusterProfiler — R でエンリッチメント解析

https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html

BiocManager::install("clusterProfiler")
library(clusterProfiler)

エンリッチメント解析には大きく二つの方法がある:

Over Representation Analysis: 何らかの解析をして得られた興味のある遺伝子群がどんな機能に偏っているかをみる。; 興味のある遺伝子/ない遺伝子が明確に二分できるとき
Gene Set Enrichment Analysis: 何らかの解析をした結果、全遺伝子に P 値などのスコアがつく場合、そのスコアの順に遺伝子を並べた上位/下位の遺伝子がどんな機能に偏っているかをみる。; 「p < 0.05 で切ってもいいけど p = 0.051 の遺伝子は本当に無関係？」みたいな問題に対処

GO enrichment analysis

GO (Gene Ontology) は、ある遺伝子がどんな機能を持っているかを共通の語彙でタグづけしたもの。もっとも大きく分けて以下の3つの分類がある:

Biological Process (BP): 遺伝子産物がどんな生物学的機能やパスウェイに属するか
Cellular Component (CC): 遺伝子産物が細胞内のどこに局在するか
Molecular Function (MF): 遺伝子産物が分子としてどういう機能をもつか

GO の情報は別途ダウンロードする。例えばヒトなら:

BiocManager::install("org.Hs.eg.db")
library(org.Hs.eg.db)

GO over-representation analysis

data(geneList, package = "DOSE")  # サンプルデータ
gene = names(geneList)[abs(geneList) > 2]

ego = enrichGO(
  gene,
  OrgDb = org.Hs.eg.db,
  keyType = "ENTREZID",
  ont = "MF",
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  universe = geneList,
  qvalueCutoff = 0.2,
  minGSSize = 10,
  maxGSSize = 500,
  readable = FALSE,
  pool = FALSE
)
ego_result = ego@result
View(ego_result)

ont: "BP", "CC", "MF", "ALL" から選ぶ。
pvalueCutoff, qvalueCutoff: それぞれ与えた P 値、Q 値以下の GO を表示する。; ここで設定してもいいけど、ここはとりあえず 1 にしておいて、 ego@result を格納してから |> dplyr::filter(qvalue < 0.1) とかする方が柔軟では。
universe: バックグラウンドの遺伝子ベクタ; 指定しないときは OrgDb の全遺伝子が使われる。
minGSSize, maxGSSize: いくつ以上/以下の遺伝子が紐づく GO までを使う。
readable: TRUE にすると ENTREZID を遺伝子シンボルに変換する。

GO Gene Set Enrichment Analysis

data(geneList, package = "DOSE")
# 渡す `geneList` はスコアの降順で並んでいる必要がある。

ego = gseGO(
  geneList,
  ont = "BP",
  OrgDb = org.Hs.eg.db,
  keyType = "ENTREZID",
  exponent = 1,
  minGSSize = 10,
  maxGSSize = 500,
  eps = 1e-10,
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  verbose = TRUE,
  seed = FALSE,
  by = "fgsea",
  scoreType = "pos"
)
ego_result = ego@result
View(ego_result)

scoreType: "pos" ならソートした上位、"neg" なら下位に濃縮する GO をみる。 (ドキュメントには載ってない？)

KEGG enrichment analysis

KEGG (Kyoto Encyclopedia of Genes and Genomes) は、主にモデル生物における遺伝子やタンパク質の分子間ネットワークに関する情報を体系化したデータベース。

利用可能な生物種を探す:

search_kegg_organism('hsa', by = 'kegg_code')
search_kegg_organism('Homo sapiens', by = 'scientific_name')

KEGG Organisms のページから探してもいい。

GO エンリッチメント解析と同様に、興味のある遺伝子/ない遺伝子で区切る解析 (over-representation) と、なんらかのスコアに基づいてソートした時に上位に濃縮するパスウェイを調べる解析 (gene set enrichment) がサポートされている。

KEGG pathway over-representation analysis

data(geneList, package = "DOSE")  # サンプルデータ
gene <- names(geneList)[abs(geneList) > 2]

kk <- enrichKEGG(
  gene,
  organism = "hsa",
  keyType = "kegg",
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  universe = geneList,
  minGSSize = 10,
  maxGSSize = 500,
  qvalueCutoff = 0.2,
  use_internal_data = FALSE
)
kk_result = kk@result
View(kk_result)

KEGG pathway gene set enrichment analysis

data(geneList, package = "DOSE")

kk <- gseKEGG(
  geneList,
  organism = "hsa",
  keyType = "kegg",
  exponent = 1,
  minGSSize = 10,
  maxGSSize = 500,
  eps = 1e-10,
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  verbose = FALSE,
  use_internal_data = FALSE,
  seed = FALSE,
  by = "fgsea",
  scoreType = "pos"
)
kk_result = kk@result
View(kk_result)

`setReadable`

enrichGO() 以外の関数には readable オプションがない。

ENTREZID のままだとどの遺伝子か分かりにくい。 setReadable は org.Hs.eg.db とかから ID と gene_symbol の対応を取得して変換する。

library(org.Hs.eg.db)
kk2 = kk |> setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")

可視化

dotplot(kk2, showCategory = 10, title = "Enriched Pathways" , split = ".sign")
# とか
cnetplot(kk2, showCategory = 5, categorySize = "pvalue")
# とか

パスウェイの図が欲しければ pathview とか:

library("pathview")
pathview(
  gene.data = geneList,
  pathway.id = "hsa04110",
  species = "hsa",
  limit = list(gene = max(abs(geneList)), cpd = 1)
)