【R】トランスクリプトーム解析のためのGO解析まとめ

Pathway Analysis

RにはGene Ontology (GO)解析やGene Set Enrichment Analysis (GSEA)を行うためのパッケージがいくつもあります．

「どのパッケージを選んだらいいのか分からない…」

という方もいらっしゃるかと思います．実際私がそうでした！

この記事ではRのGO解析を行うことのできるパッケージの中で私が実際に使ってみてお勧めしたいと感じたclusterProfilerパッケージをご紹介したいと思います．

RのGO解析用パッケージ一覧はbiopapyrusにまとまっています．

目次
Sponsored Links

GOとは？GO.dbパッケージを用いた演習

GOにまだ慣れていないのでしたらまずはGO.dbを触ってみるのもいいと思います！.dbと名前にあるパッケージはそのデータベースの情報をデータフレームのような形で保存したものです．

GO.dbに対して以下の三つの関数を通してデータベースの情報にアクセスします．

  • keys()
  • columns()
  • select()

早速GO.dbの情報をみてみましょう！

library(GO.db)

columns(GO.db)
#"DEFINITION" "GOID"       "ONTOLOGY"   "TERM"    

test_keys<-head(keys(GO.db))
#"GO:0000001" "GO:0000002" "GO:0000003" "GO:0000006" "GO:0000007" "GO:0000009"

select(GO.db,keys = test_keys ,columns = c("ONTOLOGY","TERM"),keytype = "GOID")
#        GOID ONTOLOGY                                                     TERM
#1 GO:0000001       BP                                mitochondrion inheritance
#2 GO:0000002       BP                         mitochondrial genome maintenance
#3 GO:0000003       BP                                             reproduction
#4 GO:0000006       MF    high-affinity zinc transmembrane transporter activity
#5 GO:0000007       MF low-affinity zinc ion transmembrane transporter activity
#6 GO:0000009       MF                   alpha-1,6-mannosyltransferase activity

columns(GO.db)でどんな列が含まれているのかを確認します．
“DEFINITION” 　”GOID” “ONTOLOGY” “TERM” の4つの列があります．

次にkeys()を使ってデータベースの検索を行うためのkeyを取り出しています．
GO.dbでは”GO:0000001″ のように”GO:”で始まるGOIDで検索をかけます．

最後にselect()を使ってkeyに対応するGO TERMを抽出します．

TERMはそのGOIDに対応する名前，ONTOLOGYはGO TERMの分類です．ONTOLOGYは以下の三つが存在します．
GO.dbにはBPのように二文字の略語で登録されています．

  • biological process（BP）（生物学的プロセス）
  • cellular component（CC）（細胞の構成要素）
  • molecular function（MF）（分子機能）

GOは階層構造を持っていて，上記の三つの分類の下にさらに細かい定義を持つGOタームが所属するという構成になっています．

GOの階層構造についてはbiopapyrusが詳しいです．

以上がGOの簡単な説明になります．

 

clusterProfilerによるGO解析

clusterProfilerは私が最もよく使うGO解析のパッケージです．

clusterProfilerの強みはなんと言っても解析の手軽さと豊富な可視化方法です！

Rによるパスウェイ解析の記事でもclusterProfilerをご紹介していますのでそちらもご覧ください．

clusterProfilerで行うことのできるGO解析は2種類あります．

  • GO over-representation test (GOエンリッチメント解析)
  • GO Gene Set Enrichment Analysis
Sponsored Links

GO over-representation test (GOエンリッチメント解析)

GO over-representation testは一般にGOエンリッチメント解析と呼ばれるものです．
興味のある遺伝子群（発現変動遺伝子群など）の中にどんなGOタームが濃縮されて含まれているか（エンリッチされているか）を計算する手法です．

clusterProfilerの解析演習を参考にして，DOSEパッケージのgeneListデータセットをサンプルに用いて実際に解析を行ってみましょう！

clusterProfilerではenrichGOという関数でGOエンリッチメント解析を行います．

library(clusterProfiler)
library(DOSE)
library(org.Hs.eg.db)

#geneListデータセットから発現変動遺伝子のGene IDを抽出
data(geneList)
gene <- names(geneList)[abs(geneList) > 2]

#GOエンリッチメント解析を実行
ego_result <- enrichGO(gene=gene,
                 universe=names(geneList),
                 OrgDb=org.Hs.eg.db,
                 ont="CC", 　　　　　　　　　　　　#"BP","CC","MF","ALL"から選択
                 pAdjustMethod="BH",
                 pvalueCutoff=0.01,
                 qvalueCutoff=0.05, 
                 readable=TRUE) 　　　　　　　　　　　　#Gene IDを遺伝子名に変換
 head(as.data.frame(ego_result))
#                   ID                              Description GeneRatio
#GO:0000779 GO:0000779 condensed chromosome, centromeric region    14/197
#GO:0000776 GO:0000776                              kinetochore    15/197
#GO:0000775 GO:0000775           chromosome, centromeric region    17/197
#GO:0000793 GO:0000793                     condensed chromosome    17/197
#GO:0005819 GO:0005819                                  spindle    21/197
#GO:0098687 GO:0098687                       chromosomal region    21/197

ほんの数行で解析を行うことができました！
コードの中身を見ていきましょう！

enrichGO関数は遺伝子群をEntrez Gene IDのベクトルとして受け取ります．

universeはbackground genesを定義する引数です、background genesは一般的には全遺伝子の集合を渡します.（例えばRNA-seqで検出された全遺伝子など）今回はgeneListに含まれている全ての遺伝子をbackground genesとして定義しました！

enrichGO関数の引数であるontには”BP”,”CC”,”MF”の3種類のONTOLOGYのどれを用いて解析を行うか？を渡します．readableをTRUEにするとGene IDを人にも読みやすい遺伝子名に変換してくれます．

コードの中身を見ていきましょう！

enrichGO関数は遺伝子群をEntrez Gene IDのベクトルとして受け取ります．
enrichGO関数の引数であるontには”BP”,”CC”,”MF”の3種類のONTOLOGYのどれを用いて解析を行うか？を渡します．readableをTRUEにするとGene IDを人にも読みやすい遺伝子名に変換してくれます．

このGOエンリッチメント解析の結果ですが，多くの場合似たような結果が多く含まれています．上の例だとほとんどがchromosomeに関わるものとなっています．

GOに登録されているタームにはほとんど一緒だけど少しだけ異なる機能を担うタームが数多く登録されている（Redundancyを持つ）ため，このような現象がよくみられます．
そこで，似たような結果を除外してくれるsimplifyという関数が用意されています．

ego_result.simple.simple<-simplify(ego_result)
head(as.data.frame(ego_result.simple))
#ID                              Description GeneRatio
#GO:0000779 GO:0000779 condensed chromosome, centromeric region    14/197
#GO:0000776 GO:0000776                              kinetochore    15/197
#GO:0000793 GO:0000793                     condensed chromosome    17/197
#GO:0005819 GO:0005819                                  spindle    21/197
#GO:0030496 GO:0030496                                  midbody    13/197
#GO:0031012 GO:0031012                     extracellular matrix    18/197

さきほどの結果からchromosomeに関連するものが減り，midbodyやextracellular matrixが現れました！
GO解析の結果からどんな生命現象が起こっているのかの考察を行うためには，詳細なタームをたくさん返されるより，少数の重要なタームを返してくれた方が処理がしやすいです．そのためsimplifyはとても便利な関数です！

GOエンリッチメント解析の可視化

 

#棒グラフ
barplot(ego_result.simple, drop=TRUE, showCategory=12)

 

#ドットプロット
clusterProfiler::dotplot(ego_result.simple)

#エンリッチメントマップ
clusterProfiler::emapplot(ego_result.simple)

 

#cnetプロット
clusterProfiler::cnetplot(ego_result.simple, categorySize="pvalue", foldChange=geneList)

#GO DAG graph
goplot(ego_result.simple)

 

Sponsored Links

GO Gene Set Enrichment Analysis (GSEA)

上記のGOエンリッチメント解析は興味のある遺伝子集団の機能を特定するのに非常に有用な手法です．しかし，いくつか欠点もあります．

  • 発現変動遺伝子の数が少ないデータセットからは一貫性のある生物学的な示唆を得ることが難しく，解釈が実験回ごとに，研究者ごとに異なることがあリます．
  • 細胞内プロセスは特定の遺伝子集団を協調的に変動させることが頻繁にある．例えば，特定の代謝パスウェイ中の全遺伝子が20％発現増加するとする（これは凄まじい変動と言えるでしょう）．この時に発現変動遺伝子をフォールドチェンジ2倍以上などの基準で抽出してきていると，GOエンリッチメント解析でこの代謝パスウェイを検出するのはできません．

GSEAは上記のover-representation testの欠点を補間することを目指して開発された手法です．

国立がん研のページにめちゃくちゃわかりやすいGSEAの説明が書かれていますのでぜひご覧ください！GSEAについて調べたいなら全人類これ見ればいいくらいに思っています．

ざっくり説明すると，GSEAはマイクロアレイなどの実験から得られた全ての遺伝子の2群間の発現変動量に順位をつけて一列に並べ，そのランキングのなかで特定のGOタームが発現増加側（または発現低下側）に偏っているかを検定する手法です．

GSEAもclusterProfilerによって行うことができます！

library(clusterProfiler)
library(DOSE)
library(org.Hs.eg.db)
#library(ggplot2)
data(geneList,package = "DOSE")

gse_result<- gseGO(geneList=geneList,
         OrgDb=org.Hs.eg.db,
         ont="BP",　　#"BP","CC","MF","ALL"から選択
         nPerm=1000,
         minGSSize= 120,
         pvalueCutoff=0.05,
         verbose=FALSE)

head(as.data.frame(gse_result))

 

GSEAの可視化

ridgeplot(gse_result,showCategory = 8)

gseaplot(gse_result,geneSetID = "GO:0048598")

 

参考ページ

https://bi.biopapyrus.jp/pathway/go/gene-ontology.html
https://www.ncc.go.jp/jp/ri/division/rare_cancer_research/labo/20181011150525.html
https://bioc.ism.ac.jp/packages/3.7/bioc/vignettes/enrichplot/inst/doc/enrichplot.html
https://bioinformatics-core-shared-training.github.io/cruk-summer-school-2018/RNASeq2018/html/06_Gene_set_testing.nb.html
https://bioc.ism.ac.jp/packages/3.3/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html
https://yulab-smu.github.io/clusterProfiler-book/chapter5.html

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

Sponsored Links
SKJブログはTwitterでも役に立つ情報を発信しています！

Pathway Analysisの最新記事8件