library(GSVA)
library(limma)
library(GSEABase)
geneset<-getGmt("hy_g.gmt")#可自定义，也可从数据库直接下载【Fig1】https://www.gsea-msigdb.org/gsea/msigdb/index.jsp
df<-read.table("df",header=T,row.names=1,sep="\t")#读入表达矩阵，可先进行简单筛选，例如去除表达质量差的基因
mt<-as.matrix(df)#要转为矩阵才能进行
gsvaOut<-gsva(mt,geneset)#开始计算，结果是每行一个通路，每列一个样本
#anno_col<-data.frame(sampleType=factor(c(rep("E",length(E)),rep("C",length(C)))))
#row.names(anno_col)<-colnames(df)#添加不同组别注释条
anno_col<-read.table("group",header=T,row.names=1,sep="\t")#可选，添加注释条，包含样本的组别，注意行名要和df的列名一致【Fig2】
#热图
library(pheatmap)
pheatmap(gsvaOut,annotation_col=anno_col,show_colnames=F)#【Fig3】
#箱型图
a<-gsvaOut
a$pathway<-row.names(a)
library(reshape2)
b<-melt(a,id.vars="pathway")#变为ggplot需要的长序列,按照pathway列进行转换
b$group<-anno_col[b[,2],]#添加组别
library(ggplot2)
library(ggsignif)
ggplot(b,aes(x=group,y=variable,fill=group))+geom_boxplot()+geom_point(position=position_jitterdodge())+geom_signif(comparisons=list(c("P","T")),test=wilcox.test,map_signif_level=TRUE)#comparisons参数可设置进行显著性检验的组，list格式；最后一个参数map_signif_level设为TRUE，直接以*显示，FALSE显示出具体的p值【Fig4】
Fig1-geneset
Fig1-geneset
Fig2-anno_col
Fig3-pheatmap
Fig4-boxplot
Fig4-boxplot
