GSEA

数据准备

我们需要的是带有基因id的log2foldchgange的向量，我们从de_result里面提取出来就行了

de_result$entrezid <-mapIds(org.Mm.eg.db, #.db是这个芯片数据对应的注释包
                keys=de_result$id,
                 column="ENTREZID", #clolumns参数是你要转换的ID类型是什么，只能选择一个。
                 keytype="ENSEMBL" )
#先将de_result里面的ENSEMBL转换为entrezid

geneList <- de_result$log2FoldChange
#提取所有基因的log2FoldChange
names(geneList) <- de_result$entrezid
#用对应基因的entrezid对geneList命名
geneList <- sort(geneList, decreasing = T)
#将其按照降序排列

GSEA分析

1. KEGG
   这是kegg的富集分析

require(enrichplot)
kk <- gseKEGG(geneList, organism = "mmu")

然后出现了如下报错

preparing geneSet collections...
GSEA analysis...
no term enriched under specific pvalueCutoff...
Warning message:
In fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize = minGSSize,  :
  There are duplicate gene names, fgsea may produce unexpected results

它说设置pvalueCutoff并且这里有重复的基因名，运行的话会产生意想不到的结果
于是我就手动设置pvalueCutoff

kk <- gseKEGG(geneList, organism = "mmu",pvalueCutoff = 0.05)
GSEA analysis...
leading edge analysis...
done...
Warning message:
In fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize = minGSSize, :
There are duplicate gene names, fgsea may produce unexpected results

还是报错说，有重复的基因名,在我的努力下，发现这里的问题和我是一样的
https://github.com/YuLab-SMU/clusterProfiler/issues/189,于是我按照llrs的fangfa，剔除重复的基因名

geneList <- geneList[!duplicated(names(ggeneList))]

然后再次运行

kk <- gseKEGG(geneList, organism = "mmu",pvalueCutoff = 0.05)
preparing geneSet collections...
GSEA analysis...
no term enriched under specific pvalueCutoff...

我的阈值设置得太低，导致没有结果
我们调整阈值，再来

kk <- gseKEGG(geneList, organism = "mmu",pvalueCutoff = 1)
preparing geneSet collections...
GSEA analysis...
leading edge analysis...
done...

现在就成功了！！！！
然后我们开始画图

gseKEGG_result<-as.data.frame(kk@result)
#我们先把结果提取出来，方便查看
gseaplot2(kk, "mmu04060")
#然后进行绘图

3. GO
   GO富集分析又三种,MF,BP,CC
   虽然可以用命令将他们各自求出来，但是一般情况下，我们是将他们全部求出来

kk <- gseGO(geneList,ont = "BP",OrgDb = org.Mm.eg.db)	
one of "BP", "MF", and "CC" subontologies, or "ALL" for all three
#将ont参数调整为all，就可以了
#我们可以将结果提取出来。以供查看
gseGO_result<-as.data.frame(kk@result)
View(gseGO_result)

这里我们用的是GO分析的方法，所以我们在绘图的时候是使用的GO，作为指定，而不是基因id了

gseaplot2(kk, "GO:0001816")

参考：
https://www.jianshu.com/p/a1927f75043b
https://yulab-smu.github.io/clusterProfiler-book.