一、生信分析的目的

同学们，今天来给大家分享GSEA富集分析。同样都是富集分析，它与GO和KEGG的区别在哪呢？

最主要的差别就是，GO和KEGG富集分析是针对与差异基因的分析，目的在于阐述所得到的差异基因大概会在哪些细胞通路或者功能上发挥作用。

而GSEA富集分析则是针对同学们所研究的基因进行的，它目的在于阐述目标基因可能会在哪些通路或细胞功能上发挥作用。它可以笼统的选择最相关的通路或功能展示。也可以展示自己希望展示的功能或通路，这个是最主要的功能。（那么怎么确定自己要选哪个通路展示呢，其实主要是看导师/老板的实验室有哪些蛋白，或者是自己研究的疾病与哪个通路相关性高）。

二、具体的代码

remove(list = ls()) ##清空当前环境
setwd("D:/生信代码复现/21.GSEA富集分析/21.GSEA富集分析") ##💚💚💚设置路径

#GO
#引用包
library(limma)
library(org.Hs.eg.db)
library(clusterProfiler)
library(enrichplot)

gene="CRTAC1"              #💚💚💚基因名称，我们想要研究的基因的名称
expFile="combined_RNAseq_counts.txt"     #💛💛💛表达数据文件
gmtFile="c5.go.v7.4.symbols.gmt"     #💛💛💛基因集文件

#读取文件,并对输入文件进行整理
rt=read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
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)
data=data[rowMeans(data)>0,]

#删掉正常样品
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
data=data[,group==0]
data=t(data)
rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(data))
data=t(avereps(data))

#高低表达组比较，得到logFC
dataL=data[,data[gene,]<=median(data[gene,]),drop=F]
dataH=data[,data[gene,]>median(data[gene,]),drop=F]
meanL=rowMeans(dataL)
meanH=rowMeans(dataH)
meanL[meanL<0.00001]=0.00001
meanH[meanH<0.00001]=0.00001
logFC=log2(meanH)-log2(meanL)
logFC=sort(logFC,decreasing=T)
genes=names(logFC)

#读入基因集文件
gmt=read.gmt(gmtFile)

#富集分析
kk=GSEA(logFC, TERM2GENE=gmt, pvalueCutoff = 1)
kkTab=as.data.frame(kk)
kkTab=kkTab[kkTab$pvalue<0.05,]
#kkTab=kkTab[kkTab$p.adjust<0.05,]
write.table(kkTab,file="GSEA.result-GO.txt",sep="\t",quote=F,row.names = F)
	
#输出富集的图形
termNum=5    #展示前5个通路
if(nrow(kkTab)>=termNum){
	showTerm=row.names(kkTab)[1:termNum]
	gseaplot=gseaplot2(kk, showTerm, base_size=8, title="")
	pdf(file="GSEA-GO.pdf", width=10, height=8)
	print(gseaplot)
	dev.off()
}

#输出自己想要的通路
my=read.table("my.txt", header=T, sep="\t", check.names=F)
my=as.matrix(my)
rownames(my)=my[,1]
mys=my[,2:ncol(my)]
showmy=row.names(mys)
myplot=gseaplot2(kk, showmy, base_size=8, title="")
pdf(file="GSEA-GO-myself.pdf", width=10, height=8)
print(myplot)
dev.off()





#####KEGG#####
gmtFile="c2.cp.kegg.v7.4.symbols.gmt"     #💛💛💛基因集文件
#读入基因集文件
gmt=read.gmt(gmtFile)

#富集分析
kk=GSEA(logFC, TERM2GENE=gmt, pvalueCutoff = 1)
kkTab=as.data.frame(kk)
kkTab=kkTab[kkTab$pvalue<0.05,]
#kkTab=kkTab[kkTab$p.adjust<0.05,]
write.table(kkTab,file="GSEA.result-KEGG.txt",sep="\t",quote=F,row.names = F)

#输出富集的图形
termNum=5    #展示前5个通路
if(nrow(kkTab)>=termNum){
  showTerm=row.names(kkTab)[1:termNum]
  gseaplot=gseaplot2(kk, showTerm, base_size=8, title="")
  pdf(file="GSEA-KEGG.pdf", width=10, height=8)
  print(gseaplot)
  dev.off()
}

#输出自己想要的通路
my=read.table("my.txt", header=T, sep="\t", check.names=F)   #💛💛💛输入自己想要的信号通路
my=as.matrix(my)
rownames(my)=my[,1]
mys=my[,2:ncol(my)]
showmy=row.names(mys)
myplot=gseaplot2(kk, showmy, base_size=8, title="")
pdf(file="GSEA-KEGG-myself.pdf", width=10, height=8)   #💜💜💜最终就生成了我们想要的这个GSEA富集信号通路
print(myplot)
dev.off()

三、最终的效果

①R语言中的效果：

②文件效果：

四、遇到的问题

报警：
Warning messages:
1: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, :
There are ties in the preranked stats (18.48% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.
2: In serialize(data, node载入时可能无用con) : 载入时'package:stats'可能无用
4: In serialize(data, node载入时可能无用con) : 载入时'package:stats'可能无用
6: In serialize(data, node载入时可能无用con) : 载入时'package:stats'可能无用

报错：
Error in if (abs(max.ES) > abs(min.ES)) { :
需要TRUE/FALSE值的地方不可以用缺少值

五、参考

①具体参考代码——公众号《叉叉滴同学的生信笔记》