如何使用clusterProfiler包利用eggnog-mapper軟件注釋結果做GO和KEGG富集分析

發布時間：2021-09-19 19:29:59 來源：億速云閱讀：1356 作者：小新欄目：大數據

# 如何使用clusterProfiler包利用eggnog-mapper軟件注釋結果做GO和KEGG富集分析

## 摘要
本文詳細介紹了如何利用R語言中的clusterProfiler包對eggnog-mapper軟件的注釋結果進行基因本體論(GO)和京都基因與基因組百科全書(KEGG)通路富集分析。文章包含完整的分析流程、代碼示例、結果解讀以及常見問題解決方案，適用于生物信息學研究人員進行功能基因組學研究。

## 引言

功能富集分析是解讀高通量組學數據的重要方法，能夠幫助研究者理解差異表達基因或目標基因集合的生物學功能。eggnog-mapper是一款高效的直系同源蛋白注釋工具，而clusterProfiler是R語言中最流行的富集分析工具之一。

本文將分步指導如何：

1. 準備eggnog-mapper注釋結果
2. 使用clusterProfiler進行GO富集分析
3. 進行KEGG通路富集分析
4. 可視化富集結果

## 1. 軟件安裝與環境準備

### 1.1 安裝R和必要包

```r
# 安裝clusterProfiler及相關包
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("clusterProfiler", "org.Hs.eg.db", "enrichplot", "DOSE"))

1.2 加載所需庫

library(clusterProfiler)
library(org.Hs.eg.db) # 根據物種選擇對應數據庫
library(enrichplot)
library(ggplot2)

注意：根據研究對象選擇適當的OrgDb，如小鼠用org.Mm.eg.db，大鼠用org.Rn.eg.db等。

2. eggnog-mapper結果預處理

2.1 eggnog-mapper輸出文件格式

典型輸出文件(annotations.emapper.annotations)包含多列信息，關鍵列包括：

query_name: 基因/蛋白ID
GOs: GO術語(如GO:0008150,GO:0003674)
KEGG_ko: KEGG通路標識(如ko00010,ko00230)

2.2 數據導入與清洗

# 讀取eggnog注釋結果
eggnog <- read.delim("annotations.emapper.annotations", 
                    header = TRUE, 
                    comment.char = "#",
                    stringsAsFactors = FALSE)

# 提取基因ID和GO注釋
gene_go <- eggnog[, c("query_name", "GOs")]
colnames(gene_go) <- c("GeneID", "GO")

# 去除無GO注釋的基因
gene_go <- gene_go[gene_go$GO != "", ]

2.3 構建GO注釋列表

# 分割多個GO項
go_list <- strsplit(gene_go$GO, ",")

# 創建數據框
go_annot <- data.frame(
  GeneID = rep(gene_go$GeneID, sapply(go_list, length)),
  GO = unlist(go_list)
)

# 去除重復項
go_annot <- unique(go_annot)

3. GO富集分析

3.1 準備基因列表

假設我們有一個感興趣的基因集：

interest_genes <- c("gene1", "gene2", "gene3", ...) # 替換為實際基因ID

3.2 GO富集分析

# 使用enricher函數進行自定義GO富集
ego <- enricher(
  gene = interest_genes,
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  universe = unique(go_annot$GeneID), # 背景基因集
  TERM2GENE = go_annot[, c("GO", "GeneID")]
)

# 查看簡要結果
head(ego)

3.3 使用OrgDb進行GO富集（替代方法）

如果基因ID能被OrgDb識別：

ego <- enrichGO(
  gene = interest_genes,
  OrgDb = org.Hs.eg.db,
  keyType = "SYMBOL", # 根據ID類型調整
  ont = "ALL", # 可指定"BP","MF"或"CC"
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  qvalueCutoff = 0.2
)

4. KEGG富集分析

4.1 提取KEGG注釋

# 從eggnog結果提取KEGG信息
gene_kegg <- eggnog[, c("query_name", "KEGG_ko")]
gene_kegg <- gene_kegg[gene_kegg$KEGG_ko != "", ]

# 分割多個KEGG條目
kegg_list <- strsplit(gene_kegg$KEGG_ko, ",")
kegg_annot <- data.frame(
  GeneID = rep(gene_kegg$query_name, sapply(kegg_list, length)),
  KEGG = unlist(kegg_list)
)

4.2 KEGG富集分析

ekegg <- enricher(
  gene = interest_genes,
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  universe = unique(kegg_annot$GeneID),
  TERM2GENE = kegg_annot[, c("KEGG", "GeneID")]
)

# 或使用clusterProfiler內置KEGG數據庫
ekegg <- enrichKEGG(
  gene = interest_genes,
  organism = "hsa", # 人類，其他物種需更改
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH"
)

5. 結果可視化

5.1 條形圖

barplot(ego, showCategory = 20)

5.2 點圖

dotplot(ego, showCategory = 15)

5.3 網絡圖

cnetplot(ego, categorySize = "pvalue", foldChange = geneList)

5.4 GO有向無環圖

goplot(ego)

5.5 KEGG通路圖

# 需要pathview包
library(pathview)
pathview(gene.data = geneList, pathway.id = "hsa04110")

6. 高級分析與結果解讀

6.1 GO語義相似性分析

# 簡化冗余GO項
ego2 <- simplify(ego, cutoff = 0.7, by = "p.adjust")

6.2 多組比較

# 比較不同基因集的富集結果
ck <- compareCluster(geneCluster = list(Group1 = genes1, Group2 = genes2),
                   fun = "enrichGO",
                   OrgDb = org.Hs.eg.db)
dotplot(ck)

6.3 結果導出

# 保存結果表格
write.csv(ego@result, "GO_enrichment_results.csv")
write.csv(ekegg@result, "KEGG_enrichment_results.csv")

# 保存圖形
ggsave("GO_dotplot.png", dotplot(ego), width = 10, height = 8)

7. 常見問題解答

Q1: 如何解決”No gene can be mapped”錯誤？

A: 檢查基因ID類型是否正確，必要時使用bitr()函數轉換ID：

gene_df <- bitr(interest_genes, 
               fromType = "SYMBOL", 
               toType = "ENTREZID",
               OrgDb = org.Hs.eg.db)

Q2: 如何選擇合適的背景基因集？

A: 背景集應包含所有可能被檢測到的基因，通常使用: - 整個轉錄組測序檢測到的基因 - 芯片上所有的探針對應基因 - 整個基因組基因集

Q3: 富集結果不顯著怎么辦？

A: 嘗試: 1. 放寬p值閾值(如0.1) 2. 使用更寬松的校正方法(“none”或”fdr”) 3. 檢查基因ID轉換是否正確 4. 考慮使用GSEA方法

8. 結論

本文介紹了基于eggnog-mapper注釋結果進行GO和KEGG富集分析的完整流程。通過clusterProfiler包，研究人員可以高效地解讀基因集的生物學功能和通路特征。該方法可廣泛應用于各種組學數據的生物功能分析。

參考文獻

Wu T, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021.
Huerta-Cepas J, et al. eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 2019.
Kanehisa M, et al. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017.

附錄：完整代碼示例

# 完整示例代碼
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)

# 1. 讀取數據
eggnog <- read.delim("annotations.emapper.annotations", header = TRUE, comment.char = "#")

# 2. 準備GO注釋
gene_go <- eggnog[, c("query_name", "GOs")]
gene_go <- gene_go[gene_go$GOs != "", ]
go_list <- strsplit(gene_go$GOs, ",")
go_annot <- data.frame(
  GeneID = rep(gene_go$query_name, sapply(go_list, length)),
  GO = unlist(go_list)
)

# 3. 富集分析
interest_genes <- c("gene1", "gene2", "gene3") # 替換為實際基因
ego <- enricher(gene = interest_genes,
               pvalueCutoff = 0.05,
               universe = unique(go_annot$GeneID),
               TERM2GENE = go_annot)

# 4. 可視化
dotplot(ego, showCategory = 15)

提示：實際分析時請根據具體數據調整參數和步驟。 “`

向AI問一下細節