还在为单基因GSEA发愁？看这里！果叔教你如何掌握单基因GSEA通路富集！!

library(ggplot2)#BiocManager::install("limma")library(limma)library(pheatmap)library(ggsci)lapply(c('clusterProfiler','enrichplot','patchwork'), function(x) {library(x, character.only = T)})library(org.Hs.eg.db)library(patchwork)

#输入文件expFile="after_merge.txt"sgene="GINS2"       #这里输入进行单基因GSEA的基因名称下面我们对表达矩阵处理一下，先读入rt=read.table(expFile,sep="\t",header=T,check.names=F)rt=as.matrix(rt)rownames(rt)=rt[,1]#把矩阵第一列放在R中的第一列exp=rt[,2:ncol(rt)]dimnames=list(rownames(exp),colnames(exp))data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)data=avereps(data)

#按基因表达分组group <- ifelse(data[c(sgene),]>median(data[c(sgene),]), "High", "Low")   group <- factor(group,levels = c("High","Low"))我们这里是是没有做差异分析的矩阵，我们需要找到与我们想要基因相关的，并且我们的基因也需要有差异才可以真实性接下来：#差异分析design <- model.matrix(~0+group)colnames(design) <- levels(group)fit <- lmFit(data,design)cont.matrix<-makeContrasts(High-Low,levels=design)fit2 <- contrasts.fit(fit, cont.matrix)fit2 <- eBayes(fit2)deg=topTable(fit2,adjust='fdr',number=nrow(data))Diff=deg#保存单基因分组的所有基因差异结果DIFFOUT=rbind(id=colnames(Diff),Diff)write.table(DIFFOUT,file=paste0("1.","DIFF_all.xls"),sep="\t",quote=F,col.names=F)

#展示差异最大的前30个基因Diff=Diff[order(as.numeric(as.vector(Diff$logFC))),]diffGene=as.vector(rownames(Diff))diffLength=length(diffGene)afGene=c()if(diffLength>(60)){  afGene=diffGene[c(1:30,(diffLength-30+1):diffLength)]}else{  afGene=diffGene}afExp=data[afGene,]我们对分组进行标记，#分组标签Type1=as.data.frame(group)Type1=Type1[order(Type1$group,decreasing = T),,drop=F]Type=Type1[,1]names(Type)=rownames(Type1)Type=as.data.frame(Type)#分组标签的注释颜色anncolor=list(Type=c(High="red",Low="blue"  ))
pdf(file=paste0("1.", "DIFF_heatmap.pdf"),height=7,width=6)pheatmap(afExp[,rownames(Type1)],                                                                      #热图数据         annotation=Type,                                                            #分组         color = colorRampPalette(c("blue","white","red"))(50),     #热图颜色         cluster_cols =F,                                                           #不添加列聚类树         show_colnames = F,                                                         #展示列名         scale="row",          fontsize = 10,         fontsize_row=6,         fontsize_col=8,         annotation_colors=anncolor)dev.off()
#火山图差异标准设置，这个小伙伴自行设置，并不是一定要绘制adjP=0.05aflogFC=0.5Significant=ifelse((Diff$P.Value<adjP & abs(Diff$logFC)>aflogFC), ife#开始绘制p = ggplot(Diff, aes(logFC, -log10(P.Value)))+  geom_point(aes(col=Significant),size=4)+               #size点的大小  scale_color_manual(values=c(pal_npg()(2)[2], "#838B8B", pal_npg()(1)))+  labs(title = " ")+  theme(plot.title = element_text(size = 16, hjust = 0.5, face = "bold"))+  geom_hline(aes(yintercept=-log10(adjP)), colour="gray", linetype="twodash",size=1)+  geom_vline(aes(xintercept=aflogFC), colour="gray", linetype="twodash",size=1)+  geom_vline(aes(xintercept=-aflogFC), colour="gray", linetype="twodash",size=1)#查看，不添加标记可以直接保存p#添加基因点标记，按照,可自行根据差异分析的结果进行标记point.Pvalue=0.01point.logFc=1.5其实以上果叔也教了小伙伴火山图的绘制我们直接读取结果，用下面结果去分析logFC_t=0deg$g=ifelse(deg$P.Value>0.05,'stable',             ifelse( deg$logFC > logFC_t,'UP',                     ifelse( deg$logFC < -logFC_t,'DOWN','stable') ))table(deg$g)deg$symbol=rownames(deg)df <- bitr(unique(deg$symbol), fromType = "SYMBOL",           toType = c( "ENTREZID"),           OrgDb = org.Hs.eg.db)DEG=degDEG=merge(DEG,df,by.y='SYMBOL',by.x='symbol')data_all_sort <- DEG %>%   arrange(desc(logFC))

上述就是将差异的基因绘制成火山图，而我们单独做的基因也包括在内geneList = data_all_sort$logFC #把foldchange按照从大到小提取出来names(geneList) <- data_all_sort$ENTREZID #给上面提取的foldchange加上对应上ENTREZIDhead(geneList)上述就是差异分析的过程，接下来就是我们的重头戏：#开始富集分析kk2 <- gseKEGG(geneList     = geneList,               organism     = 'hsa',               nPerm        = 10000,               minGSSize    = 10,               maxGSSize    = 200,               pvalueCutoff = 0.05,               pAdjustMethod = "none" )class(kk2)colnames(kk2@result)kegg_result <- as.data.frame(kk2)rownames(kk2@result)[head(order(kk2@result$enrichmentScore))]af=as.data.frame(kk2@result)write.table(af,file=paste0("2.","all_GSEA.xls"),sep="\t",quote=F,col.names=T)#这里我们把结果保存一下

#排序后分别取GSEA结果的前5个和后5个num=5pdf(paste0("2.","down_GSEA.pdf"),width = 8,height = 8)gseaplot2(kk2, geneSetID = rownames(kk2@result)[head(order(kk2@result$enrichmentScore),num)])dev.off()pdf(paste0("2.","up_GSEA.pdf"),width = 8,height = 8)gseaplot2(kk2, geneSetID = rownames(kk2@result)[tail(order(kk2@result$enrichmentScore),num)])dev.off()这里将上调的和下调的分析出图

也可将所以的都放在一起#排序后取前5个和后5个一起展示num=5pdf(paste0("2.","all_GSEA.pdf"),width = 10,height = 10)gseaplot2(kk2, geneSetID = rownames(kk2@result)[c(head(order(kk2@result$enrichmentScore),num),tail(order(kk2@result$enrichmentScore),num))])dev.off()

gseaplot2(kk2,          title = "Th1 and Th2 cell differentiation",  #设置标题"hsa04658", #绘制hsa04658通路的结果，通路名称与编号对应          color="red", #线条颜色          base_size = 20, #基础字体的大小          subzplots = 1:3,           pvalue_table = T) # 显示p值