2021-08-10-功能富集分析（不断分析）

功能富集分析

收到测序公司的结果以后，发现一些离谱的通路富集，于是自己想重现一下这个过程，经过一番折腾，终于摸索除了测序公司给的结果的由来。

首先我是之前都是在GSEA网站上下载相应的通路的gmt文件，[GSEA官网]（http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp#C2），在上面可以下载MSigDB所收录的基因集。下下来以后跑代码分析
比如我们先下载一个上面收录的kegg的通路集合,SYMBOL和NCBI ID都可以

image.png

##exp是我们的差异分析结果
rt1<-exp[abs(exp$log2FoldChange)>=1.0 & exp$pvalue<0.05 ,]
rt1<-arrange(rt1,-log2FoldChange)
rt1$log2FoldChange<-as.numeric(rt1$log2FoldChange)
gene=rownames(rt1)
library("org.Hs.eg.db")
entrezid <- select(org.Hs.eg.db, keys=gene, columns = 'ENTREZID', keytype = 'SYMBOL')
entrezid<-entrezid[!duplicated(entrezid$SYMBOL),]
rt1<-as.data.frame(rt1)
gene=rt1[,'log2FoldChange']
names(gene)<-entrezid$SYMBOL
gene=sort(gene,decreasing = T)
setwd("D:/SXFX/SXDT1/GSEA/")
kegmt<-read.gmt('c2.cp.kegg.v7.4.symbols.gmt')
kk=GSEA(gene,TERM2GENE=kegmt, nPerm=100,pvalueCutoff = 1)
kkTab=as.data.frame(kk)
kkTab=kkTab[kkTab$p.adjust<0.05,]
KEGG<-GSEA(gene,TERM2GENE = kegmt)

跑完上面的代码以后，没有一条通路GSEA分析的p值是有意义的，当时不是很理解，因为公司给的GO和KEGG分析都富集到了比较多的通路，我觉得不应该GSEA富集不到通路，于是我开始去跑GO和KEGG分析，然后发现公司给的结果其实就是用clusterProfiler这个R包跑的，然后我对比了下基因集，发现GSEA官网上的KEGG通路只有186条，而公司给的结果中KEGG是有三百多条，于是我就去KEGG官网上查看，发现看不大董，于是去简书上搜索，不得不说还是简书好使，很快就检索到了，发现GSEA的收集到的通路其实是相对滞后的，下面是所查阅的资料以及自己的一些整理。

KEGG通路数据库下载(人种)

参考：KEGG通路下载的几种方法 、GSEA作图

# 根据参考内容的第三种下载方法下好数据后，然后自己再制作成gmt格式的文件
gmt<-data.frame()
for (i in 1:345){
  dt<-data.frame(unlist(str_split(pathwayWithGene$hgnc_symbol[i],',')))
  colnames(dt)[1]='gene'
  dt$term=unlist(str_split(pathwayWithGene$pathway_name[i],'-'))[1]
  gmt<-rbind(gmt,dt)
}
gmt<-gmt[,c('term','gene')]

换上了新的数据库以后，果不其然，GSEA分析有结果了

这里得再Q一下clusterProfiler这个包，GO、KEGG、GSEA分析还是非常好使的，数据基本是保持最新的一个状态。
其实代码非常简单，同时GO、KEGG分析和公司的如出一辙

##KEGG分析直接用character格式的NCBI ID就好了
KEGG <- enrichKEGG(gene= entrezid$ENTREZID,organism  = 'hsa', qvalueCutoff = 0.05)
kkTab=as.data.frame(KEGG)
##GSEA分析的话，除了NCBI ID格式，还得根据logFC进行从大到小的排序
rt1<-exp[abs(exp$log2FoldChange)>=1.0 & exp$pvalue<0.05 ,]
rt1<-arrange(rt1,-log2FoldChange)
rt1$log2FoldChange<-as.numeric(rt1$log2FoldChange)
gene=rownames(rt1)
library("org.Hs.eg.db")
entrezid <- select(org.Hs.eg.db, keys=gene, columns = 'ENTREZID', keytype = 'SYMBOL')
entrezid<-entrezid[!duplicated(entrezid$SYMBOL),]
rt1<-as.data.frame(rt1)
gene=rt1[,'log2FoldChange']
names(gene)<-entrezid$ENTREZID
gene=sort(gene,decreasing = T)
#GSEA作图
kk <- gseKEGG(geneList  = gene,keyType  = 'kegg',organism = 'hsa',
  nPerm = 1000,minGSSize = 10, maxGSSize = 500,pvalueCutoff = 0.05,pAdjustMethod = "BH")
KKtable<-as.data.frame(kk)
pdf('GSEA.pdf',height = 6,width = 9)
gseaplot2(kk,geneSetID = 'hsa04140',title = kk@result$Description[56],color = 'red',
          base_size = 20,subplots = 1:2,pvalue_table = T)
dev.0ff()

image.png

同时附上KEGG和GO的barplot以及dotplot的代码

pvalueFilter=0.05                  
qvalueFilter=0.05                   
# 将 symbol 对应到 entrezid
gene=rownames(diffLab)
entrezid <- select(org.Hs.eg.db, keys=gene, columns = 'ENTREZID', keytype = 'SYMBOL')
gene=entrezid$ENTREZID
colorSel="qvalue"
if(qvalueFilter>0.05){colorSel="pvalue"}
#GO富集分析
kk=enrichGO(gene = gene,OrgDb = org.Hs.eg.db, pvalueCutoff =1, qvalueCutoff = 1, ont="all", readable =T)
GO=as.data.frame(kk)
GO=GO[(GO$pvalue<pvalueFilter & GO$qvalue<qvalueFilter),]
showNum=100
if(nrow(GO)<30){showNum=nrow(GO)}
pdf(file=paste0(name,'_GO_barplot.pdf'),width = 11,height = 7)
barplot(kk, drop = TRUE, showCategory =showNum,color = colorSel)
dev.off()
#气泡图
pdf(file=paste0(name,'_GO_bubble.pdf'),width = 11,height = 7)
dotplot(kk,showCategory = showNum, orderBy = "GeneRatio", color = colorSel)
dev.off()

#kegg
kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =1, qvalueCutoff =1)   #富集分析
KEGG=as.data.frame(kk)
KEGG$geneID=as.character(sapply(KEGG$geneID,function(x)paste(rt$SYMBOL[match(strsplit(x,"/")[[1]],as.character(rt$ENTREZID))],collapse="/")))
KEGG=KEGG[(KEGG$pvalue<pvalueFilter & KEGG$qvalue<qvalueFilter),]
#柱状图
pdf(file=paste0(name,'_KEGG_barplot.pdf'),width = 10,height = 7)
barplot(kk, drop = TRUE, showCategory = showNum, color = colorSel)
dev.off()
#气泡图
pdf(file=paste0(name,'_KEGG_bubble.pdf'),width = 10,height = 7)
dotplot(kk, showCategory = showNum, orderBy = "GeneRatio",color = colorSel)
dev.off()  

image.png image.png