clusterprofiler跑GSEA时可能存在的问题(已解决
2020-06-06 本文已影响0人
一只烟酒僧
首先说明,这个问题是在一个偶然的情况下被发现的,因此可能不具备普遍性。
事情是这样的,我在昨天(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,默认算法也稍有不同。