如果你覺得GO或KEGG不夠解釋，或許可以試試GSEA
#GSEA analysis steps
#得到差異分析結(jié)果
DESeq2_DEG = as.data.frame(read.table("路徑/DESeq2_result.txt", sep = "\t", header = T, row.names = 1))
#提取log2FoldChange這列
geneList <- as.data.frame(DESeq2_DEG$log2FoldChange)
GSEA分析需要用ENTREZID！：鸟顺！
#轉(zhuǎn)換名字族铆，調(diào)用需要的R包
library(stringr)
library(clusterProfiler)
library(org.Hs.eg.db)
geneList_tr <- bitr(row.names(DESeq2_DEG),
fromType = "SYMBOL",
toType = c("ENTREZID","ENSEMBL"),
OrgDb = org.Hs.eg.db)
#將ENTREZID與差異基因進行合并
DESeq2_DEG_withsymbol <- cbind(DESeq2_DEG, row.names(DESeq2_DEG))
names(DESeq2_DEG_withsymbol )[7]<- c("SYMBOL")
new_list <- merge(DESeq2_DEG_withsymbol, geneList_tr)
new_list <- merge(new_list, geneList_tr, by = "ENSEMBL")
geneList <- new_list$log2FoldChange
names(geneList) <- geneList_tr$ENTREZID
# 最后從大到小排序叶堆，得到一個字符串
geneList <- sort(geneList,decreasing = T)
#進行GSEA富集分析
go_result <- gseGO(geneList = geneList,
ont = "BP",
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
pvalueCutoff = 0.05,
pAdjustMethod = "BH")#p值校正方法
kegg_result <- gseKEGG(
geneList,
organism = "hsa",
keyType = "kegg",
exponent = 1,
minGSSize = 10,
maxGSSize = 500,
eps = 1e-10,
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
verbose = TRUE,
use_internal_data = FALSE,
seed = FALSE,
by = "fgsea")
head(kegg_result)
dim(kegg_result)
#按照enrichment score從高到低排序包斑，便于查看富集通路
write.table(kegg_result,"路徑/GSEA富集.txt",sep = "\t",quote = F,col.names = T,row.names = F)
#畫GSEA富集圖
library(enrichplot)
gseaplot2(
kegg_result, #gseaResult object干旧，即GSEA結(jié)果
"hsa04727",#富集的ID編號
#標題
color = "green",#GSEA線條顏色
base_size = 11,#基礎字體大小
rel_heights = c(1.5, 0.5, 1),#副圖的相對高度
subplots = 1:3, #要顯示哪些副圖 如subplots=c(1,3) #只要第一和第三個圖渠欺，subplots=1#只要第一個圖
pvalue_table = FALSE, #是否添加 pvalue table
ES_geom = "line") #running enrichment score用先還是用點ES_geom = "dot"
一個例子
如果想畫好幾條在一幅圖 可以寫
paths "hsa04510", "hsa04512", "hsa04974", "hsa05410")
gseaplot2(kk, paths)
