macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

clusterProfilerを使ってGO Enrichment Analysisを行う

2022/05/23  step5を画像に差し替え

 

先日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を追加した。

bioconductor

#conda (link)
conda activate clusterprofiler
mamba install -c bioconda bioconductor-annotationhub -y

#もしくはBiocManagerを使う
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("AnnotationHub")

 

 

実行方法

1、AnnotationHubパッケージのロード。

library(AnnotationHub)

 

2、AnnotationHubオブジェクトを作成

ah <- AnnotationHub()

これでAnnotationHubオブジェクトが使える。1分ほどかかる。

 

3、 属名や種小名で利用可能なデータベースを問い合わせる。

query(ah, "homo")

 

#データの由来
unique(ah$dataprovider)
#どのようなサンプルがあるのか
unique(ah$dataprovider)

 

利用可能なメタデータ一覧、ah$まで打ってタブキーを押す。

ah$

 

4、display() 関数を使用してブラウザでAnnotationHub オブジェクトを見ることもできる。

d <- display(ah)

立ち上がった。

上の画像でポップアップ表示されているshow new windowボタンをクリックすると、デフォルトのブラウザで開かれる(停止はその右のボタンをクリックする)

 

ブラウザでは、検索ウィンドウにタイプして調べることができる。"homo"と検索

左端の列がID。自分の所有している遺伝子IDと対応しているものを選ぶ。

 

5、Enrichment Analysisを行う生物のリソースをダウンロードする。上の流れで同定したIDを指定する。

 

6、clusterProfilerパッケージのロード。

library(clusterProfiler)

 

7、ここではsample関数でランダムに100遺伝子取ってきて、enrichGO関数(link)でGO enrichment解析を行う。ランダム抽出なのでP値カットオフは設けない。Bioconductor AnnotationData Packages(link)に登録されている生物なら、OrgDb引数をそのOrgDbデータベース名にするだけで良くて、1-6のステップはスキップできる。例えばorg.Hs.eg.db(link)にすれば、ヒトゲノムに対応する "OrgDb=org.Hs.eg.db"(マニュアル参照)。

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のカテゴリなどを指定する。

 

8、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の取り扱いなど勉強になりました。ありがとうございました。


関連