############gsekegg-20221101
suppressPackageStartupMessages({library(openxlsx) #读取.xlsx文件
library(ggplot2) #柱状图和点状图
library(stringr) #基因ID转换
library(enrichplot) #GO,KEGG,GSEA
library(clusterProfiler) #GO,KEGG,GSEA
library(GOplot) #弦图，弦表图，系统聚类图
library(DOSE)
library(ggnewscale)
library(topGO) #绘制通路网络图
library(circlize) #绘制富集分析圈图
library(ComplexHeatmap) #绘制图例
})
R.utils::setOption( "clusterProfiler.download.method",'auto' ) ##防kegg出错
setwd("/data/wanglei_lab/zhangyu/2.zy_fdy/2.TBILJN_vs_MC/tbi-20221017/6.diff/1.GSEA/gsekegg")
diff_Ast<-read.csv("/data/wanglei_lab/zhangyu/2.zy_fdy/2.TBILJN_vs_MC/tbi-20221017/6.diff/1.GSEA/diff_Ast.csv", header=TRUE,row.names= 1) #读取数据时将基因列读作行名
head(diff_Ast)#查看前面几行
dim(diff_Ast) #14955
diff_Ast$SYMBOL = rownames(diff_Ast)
library(clusterProfiler)
library(org.Mm.eg.db)
# 转成ENTREZID
df <- bitr(unique(diff_Ast$SYMBOL), fromType = "SYMBOL",toType = c( "ENTREZID"),OrgDb = org.Mm.eg.db)
merged=merge(diff_Ast,df,by='SYMBOL',all=F)
head(merged)
colnames(merged)
data_all_sort <- merged %>% arrange(desc(avg_log2FC)) #排序
GSEA_input = data_all_sort$avg_log2FC #把foldchange按照从大到小提取出来
names(GSEA_input) <- data_all_sort$ENTREZID #给上面提取的foldchange加上对应上ENTREZID
head(GSEA_input) #排序好的基因序列，而且是entrezeID的形式
#进行gseKEGG富集分析
kkgsea <- gseKEGG(geneList = GSEA_input,
organism = 'mmu', #物种
nPerm = 10000,
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 1,
pAdjustMethod = "BH" )
kkx_C2 <- setReadable(kkgsea, 'org.Mm.eg.db', 'ENTREZID') %>% as.data.frame()
write.csv(kkx_C2,file = 'diff_Ast.GSEA.result20221101-222.csv',quote = F,row.names = T)
save(GSEA_KEGG,file = 'GSEA_KEGG.Rdata')
#gseKEGG输入形式：将基因按照logFC进行从高到低排序，只需要基因列和logFC
#organism：物种,KEGG分析指定物种，物种缩写索引表详见http://www.genome.jp/kegg/catalog/org_list.html
#nPerm：permutation numbers 置换的次数
#minGSSize：通路最小基因数
#maxGSSize：通路最大基因数
#pvalueCutoff：最小p值
#pAdjustMethod：p值校正方法， 一般选择 "BH" 或 "fdr",BH较严格,fdr较温和(计算的q小些) .选BH
###作图
ridgeplot(GSEA_KEGG)
gseaplot2(GSEA_KEGG,1)
gseaplot2(GSEA_KEGG,1:30)#30是根据ridgeplot中有30个富集通路得到的
heatplot() # 获取不同富集terms间的共同基因，加上表达数据也可以看上或下调方向
cnetplot() # 获取不同富集terms间的共同基因，加上表达数据也可以看上或下调方向
upsetplot() # 不同富集terms之间的共同基因
browseKEGG(kk, "hsa04934")
install.packages("pathview")
library("pathview")
pathview(gene.data = GSEA_input,
pathway.id = "mmu05150", #上述结果中的hsa04750通路
species = "mmu",
limit = list(gene=max(abs(GSEA_input)), cpd=1),
kegg.native = TRUE) #输出的图片形式为kegg
