clusterProfiler を使ってKEGG pathway Enrichment Analysisを行う

2022/05/30 タイトル修正

2024/03/06 追記

　機能エンリッチメント解析は、生命科学におけるハイスループットなオミックスデータを解釈するために極めて重要である。この種のツールは、できるだけ多くの生物について最新のアノテーションデータベースを使用することが重要になる。このような要求に応えるために、本著者らは人気の高いBioconductorパッケージであるclusterProfiler 4.0のアップデート版をここに紹介する。このパッケージは、9年前に発表されたオリジナル版と比較して、かなり強化されている。新バージョンでは、内部でサポートされているオントロジーやパスウェイ、およびユーザーから提供されたアノテーションデータやオンラインデータベースから得られたデータに基づいて、何千もの生物における機能エンリッチメント解析のためのユニバーサルインターフェースを提供している。また、dplyr と ggplot2 パッケージを拡張し、データ操作と可視化のための整然インタフェースを提供する。その他、遺伝子セットのエンリッチメント解析や、複数の遺伝子リストのエンリッチメント結果の比較などの新機能が追加されている。clusterProfiler 4.0は、多様な生物の様々なシナリオに適用されることが期待される。

概要

clusterProfilerは、数千の生物種のコーディングおよびノンコーディングゲノミクスデータの機能特性を、最新の遺伝子アノテーションを用いて探索することをサポートします。
様々なソースからの遺伝子機能アノテーションのための普遍的なインターフェースを提供するため、様々なシナリオに適用することができます。
エンリッチメント結果へのアクセス、操作、可視化のためのインターフェースを提供し、ユーザーが効率的にデータ解釈を行えるように支援します。
複数の治療法や時間帯から得られたデータセットを一度に解析・比較することができ、異なる条件下での機能的なコンセンサスや差異を容易に明らかにすることができます。

Documentation

https://guangchuangyu.github.io/software/clusterProfiler/documentation/

Reference Manual (Version 4.4.1)

https://bioconductor.org/packages/release/bioc/manuals/clusterProfiler/man/clusterProfiler.pdf

KEGG Module Enrichment Analysis

https://guangchuangyu.github.io/2016/04/kegg-module-enrichment-analysis/

マニュアルより

KEGGモジュールのhypergeometric testとGSEAの両方がclusterProfilerでサポートされるようになりました（v３から？）。KEGGに詳しくない新規ユーザーを混乱させてしまわないように、KEGGモジュールのエンリッチメントテスト用の関数enrichMKEGGとgseMKEGGを新たに作成し、KEGGパスウェイ解析専用の関数enrichKEGGとgseKEGGは元のまま残しています。

ここではKEGG pathway Enrichment Analysisの流れだけ紹介します。

インストール

condaで仮想環境を作ってrstudioでテストした（ubuntu1８使用）。

bioconductor

Github

mamba create -n clusterprofiler python=3.9 -y
conda activate clusterprofiler
mamba install -c bioconda bioconductor-clusterprofiler -y

#今回はenrichplotも使うので導入。pathviewも導入。
mamba install -c bioconda bioconductor-pathview -y
mamba install -c bioconda bioconductor-enrichplot -y

実行方法

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

library(clusterProfiler)

２（任意）　clusterProfilerによるKEGG Enrichment Analysisは、KEGG Organismに掲載されている生物に対応している。

clusterProfilerのsearch_kegg_organism関数（link）を使ってこれを確認する。KEGG organismに登録されている3文字の識別子で検索し、scientific_nameを表示させる（３文字の識別子はKEGG Organismsで確認できる）。ここでは識別子lboを指定している。

　search_kegg_organism('lbo', by='kegg_code')

出力

３（任意）、scientific_nameが分かったら、search_kegg_organism関数（link）を使ってclusterProfilerで利用可能なアノテーションを確認する。ステップ２で確認したscientific_nameを指定している。

lbo <- search_kegg_organism('Leptolyngbya boryana', by='scientific_name')
head(lbo)

ここでは１つしかヒットしないが、１つの種から複数の株や品種のアノテーションがKEGGに登録されている場合、一覧がリスト表示される（*１）。

４、#DESeq2の検定結果から始める(クラシカルな２群間比較の手順)。ORAでは、データを整形して遺伝子IDだけ取り出す。GSEA（preranked）では遺伝子IDとlog fold changeを取り出す。

dds <- DESeqDataSetFromMatrix(countData = count, colData = coldata, design = ~ condition)
dds <- DESeq(dds)
res <- results(dds)

head(res) #確認
sort <- as.data.frame(res[order(res$padj),]) #FDR順にソート、データフレームとして返す
sig <- sort[sort$padj < 0.01,] #ここではFDR0.01以下を返す
dim(sig) #遺伝子数をチェック
induce <- sig[sig$log2FoldChange > 0,] #誘導遺伝子を返す
repress <- sig[sig$log2FoldChange < 0,] #抑制遺伝子を返す
inducedgene <- rownames(induce) #誘導された遺伝子名リストを返す
head(inducedgene)
repressedgene <- rownames(repress) #誘導された遺伝子名リストを返す
head(repressedgene)

(欠損値がある場合、除去にはna.omit関数などが使える（参考）)

５、enrichKEGG関数（link）を使ってKEGG pathway over-representation analysis (ORA)を行う。入力遺伝子リスト、3文字の識別子、P値のカットオフを指定する。オプションとして、多重検定補正の方法とq値（FDR）のカットオフ値、各カテゴリのリストの下限と上限値などを指定できる。認識できる遺伝子IDは、kegg geneID以外に、ncbi-geneid、ncib-proteinid、uniprotがある（keyTypeオプションで指定）。生物により対応しているかは異なる。

#データの読み込み(KEGG geneIDのリスト)
list <- read.table("geneID.txt",sep = "\t", header = T)

#enrichKEGGによるKEGG pathway over-representation analysis。ここでは4のinducedgeneを使用
orainduce <- enrichKEGG(gene = inducedgene, organism = 'lbo', pvalueCutoff = 0.05, minGSSize = 10, maxGSSize = 500, pAdjustMethod = "BH", qvalueCutoff = 0.05)
head(orainduce)

gene a vector of entrez gene id.
organism supported organism listed in 'http://www.genome.jp/kegg/catalog/org_list.html'
keyType one of "kegg", 'ncbi-geneid', 'ncib-proteinid' and 'uniprot'
pvalueCutoff Cutoff value of pvalue.
pAdjustMethod one of "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"
universe background genes
minGSSize minimal size of genes annotated by Ontology term for testing.
maxGSSize maximal size of genes annotated for testing
qvalueCutoff qvalue cutoff
use_internal_data logical, use KEGG.db or latest online KEGG data

注；Ensembl IDを使っていてkeggが認識できない場合、ヒトなどならDAVIDで変換する（例　Ensembl gene ID => ENTREZ gene ID）。keyType = "kegg"は、 "kegg", ’ncbi-geneid’, ’ncib-proteinid’, ’uniprot’から選べるが、バージョンによっては該当するURLにアクセスできずエラーになる。

６、結果を視覚化する。ここではstep5のORAの結果を想定している。

library(enrichplot)

#barplotで視覚化(*2)
barplot(orainduce, showCategory=10) #横軸は検出語彙数
mutate(orainduce, qscore = -log(p.adjust, base=10)) %>% 
    barplot(x="qscore") #横軸はq値の-log

#dot plotで視覚化
dotplot(orainduce, showCategory=10) + ggtitle("dotplot for kegg ORA")

#upset plotで視覚化(異なる遺伝子セット間の重複を計算)
upsetplot(orainduce)

#特定のパスウェイでのアサイン率（赤色）を見る。3文字の識別子＋KEGG pathwayID(5桁の数値の識別子)で指定する
#browseKEGG関数（link）
browseKEGG(orainduce, 'lbo00195')

７、gseKEGG関数（link）を使ってKEGG pathway gene set enrichment analysis（GSEA pre-ranked）を行う（*1, 3）。

#DESeq2の結果のresオブジェクトを使う(step4の最初)
data <- as.data.frame(res) #データフレームとして返す
geneList <- data$log2FoldChange #倍率変化列をnumericベクトルとして返す
head(geneList)
names(geneList) <- as.character(rownames(data)) #名前付きベクトルとして返す
head(geneList)
geneList = sort(geneList, decreasing = TRUE) #pre-ranked GSEAを行うには降順に降順にソートされている必要がある。
head(geneList)

#gseKEGGを使ってpre-ranked GSEAを実行
kk2 <- gseKEGG(geneList = geneList, organism = 'lbo', minGSSize = 10, pvalueCutoff = 0.05, verbose = FALSE)
head(kk2)

#ドットプロットで視覚化 
dotplot(kk2, showCategory=10) + ggtitle("dotplot for kegg GSEA")

#ridgeplot(link)
ridgeplot(kk2)

#gseaplot2(link)、例えばランクの2番目をチェックする
gseaplot2(kk2, geneSetID = 2, title = kk2$Description[2])
=> プロットの上部は、遺伝子セット（例：muscle system process）のランニングESを示し、解析がランク付けされたリストを下っていく様子を示している。プロットの中央部分は、遺伝子セットのメンバーが遺伝子のランク付けされたリストのどこに表示されるかを示している。プロットの下部は、ランク付けされた遺伝子のリストを下るにつれて、ランキングメトリクスの値（例えば、フォールドチェンジ）を示している。ランキングの値は、ランク付けされたリストの下に行くに従って、正から負になる

geneList order ranked geneList
organism supported organism listed in 'http://www.genome.jp/kegg/catalog/org_list.html'
keyType one of "kegg", 'ncbi-geneid', 'ncib-proteinid' and 'uniprot'
exponent weight of each step
nPerm permutation numbers
minGSSize minimal size of each geneSet for analyzing
maxGSSize maximal size of genes annotated for testing
pvalueCutoff pvalue Cutoff
pAdjustMethod pvalue adjustment method
verbose print message or not
use_internal_data logical, use KEGG.db or latest online KEGG data
seed logical

引用

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