首先说明，这个问题是在一个偶然的情况下被发现的，因此可能不具备普遍性。
事情是这样的，我在昨天（6月4日）用clusterprofiler跑了以下语句：GSEA_GO_res<-gseGO(gene_fc,OrgDb = org.Mm.eg.db,pvalueCutoff = 1)，之后只保存了csv文件，然后关机，在今天（6月5日），由于之前忘记保存rdata，因此重新运行了原来的代码，发现在相同的代码下，两次运行的结果很不同，如下图：
这是第一次跑的，保存成csv
这是第二次跑的
当然，我也怀疑过自己是否犯了一些低级错误，如将基因集搞混之类的，因此我将两次的结果做了比较，代码如下
#这里的result_check是导入的第一次跑的结果
#第二次跑的结果已经在内存中，变量名为GSEA_GO_res
result_check<-read.csv("../chongdanyang/file/富集分析/GSEA/P5KO.VS.P5WT_GSEA_GO_res.csv")
result_check[1:6,1:5]
#我想检查一下两次跑的结果中，富集到的GOterm是否相似
check_res=c()
for (i in 1:length(GSEA_GO_res@result$Description)) {
check_res1=grepl(paste("^",as.character(GSEA_GO_res@result$Description[i]),"$",sep = ""),result_check$Description)
check_res1=ifelse(setequal(which(check_res1==TRUE),integer(0)),FALSE,TRUE)
check_res=c(check_res,check_res1)
}
table(check_res)
#check_res
#FALSE TRUE
# 12 5750
dim(result_check)
#[1] 5762 11
dim(GSEA_GO_res@result)
#[1] 5762 11
length(GSEA_GO_res@geneSets)
#16070
#数据集中总的geneset数目为16070
结果如上，两次结果大部分term（5750/5762）都是相同的，只有极少数term不同，同时发现，虽然term可以被富集到，但是两次计算的NES以及p值都不同。如nucleotide biosynthetic process这个term，在第一次的结果中，NES以及p值分别为 1.57113781 0.001148106，而在第二次的结果中，上述值为 1.569829 0.001156069。
这个故事告诉我们：跑程序一定要把所有中间文件备份好，这样你就不会发现这些糟心的问题了（笑）
知道为什么的小伙伴欢迎留言！在此先谢过！
一个心明眼亮的人的答案https://www.jianshu.com/p/2d0474f5c6f9
贴上余老师的邮件回复：用nPerm=10000，p值稳定性会比较好。不过如果你用最新版本的clusterProfiler，默认算法也稍有不同。
网友评论