首先說明，這個問題是在一個偶然的情況下被發(fā)現(xiàn)的瑟由，因此可能不具備普遍性潮瓶。
事情是這樣的，我在昨天（6月4日）用clusterprofiler跑了以下語句：GSEA_GO_res<-gseGO(gene_fc,OrgDb = org.Mm.eg.db,pvalueCutoff = 1)俱诸，之后只保存了csv文件，然后關機赊舶，在今天（6月5日），由于之前忘記保存rdata赶诊，因此重新運行了原來的代碼笼平，發(fā)現(xiàn)在相同的代碼下，兩次運行的結果很不同舔痪，如下圖：
當然，我也懷疑過自己是否犯了一些低級錯誤锄码，如將基因集搞混之類的夺英，因此我將兩次的結果做了比較晌涕，代碼如下
#這里的result_check是導入的第一次跑的結果
#第二次跑的結果已經(jīng)在內存中，變量名為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
#數(shù)據(jù)集中總的geneset數(shù)目為16070
結果如上余黎，兩次結果大部分term（5750/5762）都是相同的，只有極少數(shù)term不同载萌，同時發(fā)現(xiàn)惧财，雖然term可以被富集到，但是兩次計算的NES以及p值都不同扭仁。如nucleotide biosynthetic process這個term垮衷，在第一次的結果中，NES以及p值分別為 1.57113781 0.001148106乖坠，而在第二次的結果中搀突，上述值為 1.569829 0.001156069。
這個故事告訴我們：跑程序一定要把所有中間文件備份好熊泵，這樣你就不會發(fā)現(xiàn)這些糟心的問題了（笑）
知道為什么的小伙伴歡迎留言仰迁！在此先謝過！
一個心明眼亮的人的答案http://www.reibang.com/p/2d0474f5c6f9
貼上余老師的郵件回復：用nPerm=10000戈次，p值穩(wěn)定性會比較好轩勘。不過如果你用最新版本的clusterProfiler，默認算法也稍有不同怯邪。
