Y叔蕉绿:仅4.1%富集分析研究明确报告背景选择方式

# ==============================================# 代码功能：对比不同背景基因集对GO/KEGG富集分析结果的影响# 核心：验证全基因组/表达基因集 作为universe时的p值相关性# ==============================================rm(list = ls())# 设置R下载超时时间：将超时阈值设为极大值，彻底避免下载包/数据时因超时而中断options(timeout = 23333333333333333333333333333333333)# 加载依赖包library(ggplot2)library(clusterProfiler)  # 富集分析主包library(org.Hs.eg.db)     # 人类基因注释数据库# ---------------------- 加载示例数据 ----------------------# 获取人类全基因组所有基因（作为背景基因集）allgenes <- keys(org.Hs.eg.db)# 加载DOSE包内置的演示基因列表data(geneList, package="DOSE")# 提取前200个基因作为差异基因集（待富集基因）de <- names(geneList)[1:200]# ---------------------- 功能富集分析 ----------------------# 开启参数：强制使用自定义的背景基因集options(enrichment_force_universe=TRUE)# GO生物学过程(BP)富集分析bp_allgene <- enrichGO(de, OrgDb = 'org.Hs.eg.db', ont="BP", universe=allgenes)bp_bg1 <- enrichGO(de, OrgDb = 'org.Hs.eg.db', ont="BP", universe = names(geneList))# KEGG通路富集分析bp_allgene_kegg <- enrichKEGG(de, universe=allgenes)bp_bg1_kegg <- enrichKEGG(de, universe = names(geneList))# 关闭参数：自动将背景基因集与注释基因集取交集options(enrichment_force_universe=FALSE)# 再次执行GO/KEGG富集（自动取交集模式）bp_allgene2 <- enrichGO(de, OrgDb = 'org.Hs.eg.db', ont="BP", universe=allgenes)bp_bg2 <- enrichGO(de, OrgDb = 'org.Hs.eg.db', ont="BP", universe = names(geneList))bp_allgene2_kegg <- enrichKEGG(de, universe = allgenes)bp_bg2_kegg <- enrichKEGG(de, universe = names(geneList))# ---------------------- 可视化绘图 ----------------------# 获取富集分析专用配色方案cols <- enrichplot:::get_enrichplot_color()# 自定义绘图函数：绘制两组富集结果的p值相关性散点图# res1：第一组富集结果  res2：第二组富集结果myplot <- function(res1, res2) {# 按ID合并两组富集结果  x <- merge(res1@result, res2@result, by="ID")# 绘制散点图+拟合线+参考对角线  ggplot(x, aes(pvalue.x, pvalue.y)) +     geom_point(shape=21, fill=cols[2], size=2.5, alpha=.8, stroke=.02, color='white') +     geom_smooth(method='lm', color=cols[1]) +    geom_abline(slope=1, intercept=0, linetype='dashed', linewidth=1) +    xlim(0, 1) + ylim(0, 1) + DOSE::theme_dose()}# 绘制6张子图（GO富集3张 + KEGG富集3张）# 图1：GO-BP | 全基因组背景 vs 表达基因背景（强制模式）p1 <- myplot(bp_allgene, bp_bg1) +   xlab('All genes in the genome as universe') +  ylab("All expressed genes as universe")# 图2：GO-BP | 表达基因背景(交集) vs 表达基因背景(强制)p2 <- myplot(bp_bg2, bp_bg1) +  xlab("All expressed genes as universe\n(intersected with gene sets)") +  ylab("All expressed genes as universe")# 图3：GO-BP | 全基因组背景(交集) vs 表达基因背景(交集)p3 <- myplot(bp_allgene2, bp_bg2) +   xlab('All genes in the genome as universe\n(intersected with gene sets)') +  ylab("All expressed genes as universe\n(intersected with gene sets)")# 图4：KEGG | 全基因组背景 vs 表达基因背景（强制模式）p4 <- myplot(bp_allgene_kegg, bp_bg1_kegg) +   xlab('All genes in the genome as universe') +  ylab("All expressed genes as universe")# 图5：KEGG | 表达基因背景(交集) vs 表达基因背景(强制)p5 <- myplot(bp_bg2_kegg, bp_bg1_kegg) +  xlab("All expressed genes as universe\n(intersected with gene sets)") +  ylab("All expressed genes as universe")# 图6：KEGG | 全基因组背景(交集) vs 表达基因背景(交集)p6 <- myplot(bp_allgene2_kegg, bp_bg2_kegg) +   xlab('All genes in the genome as universe\n(intersected with gene sets)') +  ylab("All expressed genes as universe\n(intersected with gene sets)")# 组合所有子图：3列布局，自动添加A-F标签，并赋值给组合图对象combined_plot <- aplot::plot_list(p1, p2, p3, p4, p5, p6, ncol=3, tag_levels='A')# 导出组合图为高清PDF文件（适配期刊投稿标准）ggsave(  filename = "clusterProfiler_p1_6.pdf",  # 输出的PDF文件名  plot = combined_plot,                    # 指定待保存的图形  width = 14,                              # PDF宽度（英寸）  height = 8,                             # PDF高度（英寸）  device = "pdf"# 格式设置为PDF)