2022/05/23 step5を画像に差し替え
2024 03/05 dockerイメージ例追記
先日clusterProfilerを使ってKEGG termのエンリッチメント解析を行う例を紹介しました。今回はclusterProfilerを使ってGO Enrichment Analysisを行う流れを紹介します。Bioconductor AnnotationData Packages（link）として公開されているアノテーションであればそのパッケージをロードするだけで良いのですが、そうじゃない生物がほとんどなので、ここでは、annotationhubからGO termを含むアノテーションを取ってきてGO Enrichment Analysisを行います。clusterProfilerのマニュアルでも推奨されている方法になります。
GO analysis using clusterProfiler
https://guangchuangyu.github.io/2016/01/go-analysis-using-clusterprofiler/
AnnotationHub Dcocumentation
https://bioconductor.org/packages/3.14/bioc/vignettes/AnnotationHub/inst/doc/AnnotationHub.html
マニュアルより
clusterProfilerは、Gene Ontologyの過剰発現テスト（ORA）と遺伝子セットエンリッチメント解析（GSEA）をサポートしています。OrgDbオブジェクト、GMTファイル、ユーザー自身のデータからのGOアノテーションをサポートしています。github版のclusterProfilerでは、enrichGOとgseGO関数からorganismのパラメータが削除され、OrgDbのパラメータが追加されたので、OrgDbオブジェクトがある種であれば、clusterProfilerで解析ができるようになりました。Bioconductorは既に約20種のOrgDbを提供しており、http://bioconductor.org/packages/release/BiocViews.html#___OrgDb。ユーザはAnnotationHub経由でOrgDbを構築することができます。
インストール
前回condaを使って作った環境にannotationhubを追加した。
#conda (link)
conda activate clusterprofiler
mamba install -c bioconda bioconductor-annotationhub -y
#もしくはBiocManagerを使う
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("AnnotationHub")
#追記 henryc101さんのdocker image (dockerhub)(注; AnnotationHubは入ってない)
docker pull henryc101/clusterprofiler:4.10.0
#run
docker run -itv $PWD:/data -w /data --rm henryc101/clusterprofiler:4.10.0 R
実行方法
clusterProfilerの説明は6から。1-5はAnnotationHubによるアノテーションの取得の長さ。
１、AnnotationHubパッケージのロード。
library(AnnotationHub)
２、AnnotationHubオブジェクトを作成
ah <- AnnotationHub()
これでAnnotationHubオブジェクトが使える。1分ほどかかる。
３、 属名や種小名で利用可能なデータベースを問い合わせる。
query(ah, "homo")
#データの由来
unique(ah$dataprovider)
#どのようなサンプルがあるのか
unique(ah$dataprovider)
利用可能なメタデータ一覧、ah$まで打ってタブキーを押す。
ah$
４、display() 関数を使用してブラウザでAnnotationHub オブジェクトを見ることもできる。
d <- display(ah)
立ち上がった。
上の画像でポップアップ表示されているshow new windowボタンをクリックすると、デフォルトのブラウザで開かれる（停止はその右のボタンをクリックする）
ブラウザでは、検索ウィンドウにタイプして調べることができる。"homo"と検索
左端の列がID。自分の所有している遺伝子IDと対応しているものを選ぶ。
５、Enrichment Analysisを行う生物のリソースをダウンロードする。上の流れで同定したIDを指定する。
６、clusterProfilerパッケージのロード。
library(clusterProfiler)
７、ここではsample関数でランダムに100遺伝子取ってきて、enrichGO関数（link）でGO enrichment解析を行う。ランダム抽出なのでP値カットオフは設けない。Bioconductor AnnotationData Packages（link）に登録されている生物なら、OrgDb引数をそのOrgDbデータベース名にするだけで良くて、１－６のステップはスキップできる。ヒトゲノムなら "OrgDb=org.Hs.eg.db"にする（マニュアル参照）。org.Hs.eg.db（link）を指定している。
sample_gene <- sample(keys(rice), 100)
str(sample_gene) #自分のデータ解析に向けて遺伝子IDをチェックしておく
sample_test <- enrichGO(sample_gene, OrgDb=rice, pvalueCutoff=1, qvalueCutoff=1)
head(summary(sample_test))
#bar plot
barplot(sample_test)
- gene a vector of entrez gene id.
- OrgDb OrgDb
- keytype keytype of input gene
- ont One of "MF", "BP", and "CC" subontologies. #allもあり
- pvalueCutoff Cutoff value of pvalue.\
- pAdjustMethod one of "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"
- universe background genes
- qvalueCutoff qvalue cutoff
- minGSSize minimal size of genes annotated by Ontology term for testing.
- maxGSSize maximal size of genes annotated for testing
- readable whether mapping gene ID to gene Name
> head(summary(sample_test))
実際の解析では遺伝子キーを取ってきて行う。超幾何検定のP値カットオフ、多重検定の補正方法とq値カットオフ、最大サイズ、最小サイズ、GO termのカテゴリなどを指定する。
８、gseGO関数（link）を使って遺伝子セットエンリッチメント解析（GSEA）を行う。RNA seqのデータを使うなら、前回説明したように遺伝子キーとlog fold changeのオブジェクトを指定する。
gsecc <- gseGO(geneList=geneList, ont="CC", OrgDb=org.Hs.eg.db, verbose=F)
head(summary(gsecc))
引用
clusterProfiler 4.0: A universal enrichment tool for interpreting omics data
Tianzhi Wu, Erqiang Hu, Shuangbin Xu, Meijun Chen, Pingfan Guo, Zehan Dai, Tingze Feng, Lang Zhou, Wenli Tang, Li Zhan, Xiaocong Fu, Shanshan Liu, Xiaochen Bo, Guangchuang Yu
Innovation (Camb). 2021 Jul 1;2(3):100141
参考
もじゃもじゃのほうさんのQiita記事を参考にしました。IDの取り扱いなど勉強になりました。ありがとうございました。
関連
