我们在使用clusterProfiler做GO富集分析的时候，需要在Bioconductor上下载目标物种的注释包，比如人类的：

有一些物种的注释包是能够下载到的，Bioconductor上提供了以下19个物种的Org类型的包，包含了这些模式物种的GO及KEGG注释信息：

然而水稻的注释包在Bioconductor是下载不到的，有两种解决办法：

一种是利用AnnotationHub在线检索抓取OrgDb, 但是这些包是用ENTREZID，需要先将RAP-DB或者MSU转为ENTREZID才行

library(AnnotationHub)


hub <- AnnotationHub()

# 检索水稻相关 OrgDb 
query(hub, c("Oryza sativa", "OrgDb"))

# AnnotationHub with 3 records
# # snapshotDate(): 2025-10-29
# # $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
# # $species: Oryza sativa_subsp._japonica, Oryza sativa_Japonica_Group, Oryza...
# # $rdataclass: OrgDb
# # additional mcols(): taxonomyid, genome, description,
# #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
# #   rdatapath, sourceurl, sourcetype
# # retrieve records with, e.g., 'object[["AH119907"]]'
# 
# title
# AH119907 | org.Oryza_sativa_(japonica_cultivar-group).eg.sqlite
# AH119908 | org.Oryza_sativa_Japonica_Group.eg.sqlite
# AH119909 | org.Oryza_sativa_subsp._japonica.eg.sqlite

# 假设查询到水稻对应的编号为 AH119907（实际编号以你运行检索时的输出为准）
options(AnnotationHub.url = "https://mirrors.tuna.tsinghua.edu.cn/bioconductor")  #镜像加速一下
rice_orgdb <- hub[['AH119907']]
#rice_orgdb1 <- hub[['AH119909']]
#rice_orgdb2 <- hub[['AH119908']]

# 此时你的 rice_orgdb 已经加载好了
# 确认支持的ID类型（通常包含 ENTREZID, GO, EVIDENCE 等）
keytypes(rice_orgdb)


# 1. 首先获取数据库中前 6 个基因的 ENTREZID 主键
top_gene_keys <- head(keys(rice_orgdb, keytype = "ENTREZID"))

# 2. 使用 select 函数提取这 6 个基因对应的多维度注释信息
top_annotations <-  AnnotationDbi::select(
  x = rice_orgdb, 
  keys = top_gene_keys, 
  keytype = "ENTREZID", 
  columns = keytypes(rice_orgdb)
)

# 3. 在控制台展示这前几行的对照表内容
print(top_annotations)



## 获得富集分析

library(clusterProfiler)
library(dplyr)       # 用于数据清洗和表格操作

# 模拟你的 SYMBOL 基因列表（请替换为你自己的差异基因名）
set.seed(2026)
my_large_loc_genes <- sample(top_annotations$SYMBOL,1000)
# 确保去掉基因名两端的空格
my_large_loc_genes <- trimws(my_large_loc_genes)




# 运行 GO 富集分析
go_result <- enrichGO(
  gene         = my_large_loc_genes,       # 差异基因列表（必须是字符型向量）
  OrgDb        = rice_orgdb,     # 你刚刚下载的水稻 OrgDb 对象
  keyType      = "SYMBOL",     # 输入基因的 ID 类型
  ont          = "ALL",          # 可选 "BP", "CC", "MF" 或 "ALL"
  pAdjustMethod = "BH",          # 多重假设检验矫正方法
  pvalueCutoff = 0.05,           # p值显著性阈值
  qvalueCutoff = 0.2             # q值阈值
)

# 1. 查看富集分析结果表格
go_result_df <- as.data.frame(go_result)
head(go_result_df)

# 2. 可视化：柱状图与气泡图
barplot(go_result, showCategory = 20, title = "GO Enrichment - Rice")
dotplot(go_result, showCategory = 20, title = "GO Enrichment - Rice")



# 运行 KEGG 富集分析（直接在线请求 KEGG 官网最新数据）

gene_annotations <- bitr(
  geneID   = my_large_loc_genes, 
  fromType = "SYMBOL",         # 输入是 SYMBOL
  toType   = c("ENTREZID", "GO", "EVIDENCE"), # 想转换出的类型
  OrgDb    = rice_orgdb        # 你下载的水稻 rice_orgdb 对象
)

head(gene_annotations)

kegg_result <- enrichKEGG(
  gene         = gene_annotations$ENTREZID,       # 必须是 NCBI ENTREZID（纯数字）
  organism     = "osa",          # 水稻日本晴的 KEGG 缩写是 osa
  keyType      = "kegg",         # 对应纯数字 id
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH"
)

# 查看 KEGG 结果
kegg_result_df <- as.data.frame(kegg_result)
head(kegg_result_df)

# 气泡图可视化
dotplot(kegg_result, showCategory = 20, title = "KEGG Enrichment - Rice")