GO、KEGG（批量分组）分析及可视化

这篇帖子其实是更新、补充、解决一下问题的。我们号写过很多GO、KEGG富集分析的可视化，数不胜数，可以在公众号检索“富集”了解更多。我们演示的时候都是直接提供了富集的结果文件，一般演示为了图方便，也是利用在线工具cytoscape做的。结果一伙伴最近提问有做过GO、KEGG富集的R语言帖子么，突然发现这样的内容还没有正经写过，所以这里补充一下。

1、GO、KEGG分析

首先我们做一下单独的GO、KEGG分析，这里我们使用的是引用很高的，基本上人人都在用的余老师的R包-clusterProfiler，相信大家都很熟悉了。我们这里创新在于使用循环，对上下调基因分别进行GO、KEGG分析，并保存结果。初学者小伙伴可以好好理解下代码意思。

setwd("D:/KS科研分享与服务公众号/GO和KEGG分析")
library(clusterProfiler)
library(ggplot2)
library(org.Hs.eg.db)
#筛选显著性上下调基因
df <- read.csv('DEGs_trans.csv', header = T)
df_up <- df[which(df$log2FoldChange>2 & df$padj < 0.01),]#至于什么阈值自己定
df_down <- df[which(df$log2FoldChange< -2 & df$padj < 0.01),]
#做GO和KEGG分析需要进展gene id转化

sig_gene <- list(df_up, df_down)
names(sig_gene) <- c("df_up","df_down")

#GO和KEGG结果存储
sig_gene_GO <- list()
sig_gene_KEGG <- list()

#里面分析阈值请自己确定
for (i in 1:length(sig_gene)){
  
  entrezid_all = mapIds(x = org.Hs.eg.db,  #按照自己的物种修改
                        keys = sig_gene[[i]]$gene, 
                        keytype = "SYMBOL", 
                        column = "ENTREZID")
  entrezid_all  = na.omit(entrezid_all) 
  entrezid_all = data.frame(entrezid_all) 
  
  GO_enrich = enrichGO(gene = entrezid_all[,1],
                       OrgDb = org.Hs.eg.db, #按照自己的物种修改
                       keyType = "ENTREZID", 
                       ont = "BP", #可以选择All、或者CC、MF
                       pAdjustMethod = "BH",
                       pvalueCutoff = 0.05,
                       qvalueCutoff = 0.05,
                       readable = T) 
  GO_enrich  = GO_enrich@result
  
  sig_gene_GO[[i]] <- GO_enrich
  names(sig_gene_GO)[i] <- names(sig_gene)[i]
  
  
  KEGG_enrich = enrichKEGG(gene = entrezid_all[,1], 
                           keyType = "kegg",
                           pAdjustMethod = 'BH',  
                           organism= "hsa",  #https://www.kegg.jp/brite/br08611 #按照自己的物种修改
                           qvalueCutoff = 0.05,
                           pvalueCutoff=0.05)
  KEGG_enrich  = KEGG_enrich@result
  
  sig_gene_KEGG[[i]] <- KEGG_enrich
  names(sig_gene_KEGG)[i] <- names(sig_gene)[i]
  
}


#保存为csv文件
for (i in 1:2){
  A <- sig_gene_GO[[i]]
  write.csv(A,file =paste(names(sig_gene)[i],"_GO.csv"))
  
}

for (i in 1:2){
  A <- sig_gene_KEGG[[i]]
  write.csv(A,file =paste(names(sig_gene)[i],"_KEGG.csv"))
  
}

最后可视化我们决定不再使用之前的哪些可视化方式了，因为写过太多了，再写就没意思了。正好最近看到老俊俊的生信笔记分享的一个富集可视化R包aPEAR，我们就利用这个包顺便做一下可视化。是网络的形式，GO、KEGG结果都可以展示，还是可以。更多详细的用法请查看这个包的帮助文档！

install.packages("aPEAR")#https://github.com/ievaKer/aPEAR
library(aPEAR)
#相关阈值设定清查阅读 ?enrichmentNetwork
#enrichmentNetwork返回的是ggplot对象，所以可以用ggplot2主题修饰
enrichmentNetwork(sig_gene_KEGG[[1]], 
                  colorBy = 'pvalue', 
                  colorType = 'pval')+
  scale_color_gradientn(colours = c("#B83D3D",'white','#1A5592'),
                        name = "pvalue")

image.png

到这里我们这篇帖子还没结束，有两个原因。第一我们之前总是说一个函数，多组GO分析（例如：Clusterprofile多分组富集分析及可视化），但是没有说过KEGG，有小伙伴去做的时候出错，其实参数里面选择KEGG，设置对就可以进行这个分析。第二个问题是有小伙伴发来图让复现，是富集结果的展示，乍一看很复杂，既是网络图，又是多组的，其实很简单，clusterProfiler多组富集分析和enrichplot早就解决了这个问题。所以这两个问题我们归为一个进行解决。

image.png

（reference：A single-cell transcriptomic atlas of exercise-induced antiinflammatory and geroprotective effects across the body）

2、分组GO、KEGG分析

我们还是利用之前的分组基因的文件，进行多组KEGG分析（这里就不再展示GO)，之后进行可视化。首先是分析，很简单：

setwd("D:/KS科研分享与服务公众号/多组富集分析")
library(Seurat)
library(SeuratData)
library(clusterProfiler)
library(enrichplot)
df_sig <- read.csv("df.csv", header = T)
#数据如此格式即可，其他的数据整理成此格式即可
group <- data.frame(gene=df_sig$gene,group=df_sig$cluster)#分组情况

#gene转化为ID
Gene_ID <- bitr(df_sig$gene, fromType="SYMBOL", 
                toType="ENTREZID", 
                OrgDb="org.Hs.eg.db")

#构建文件并分析
data  <- merge(Gene_ID,group,by.x='SYMBOL',by.y='gene')
data_KEGG <- compareCluster(ENTREZID~group,
                            data=data,
                            fun = "enrichKEGG",#函数选择什么定义什么分析
                            pAdjustMethod = "BH",
                            pvalueCutoff = 0.05,
                            qvalueCutoff = 0.05,
                            organism= "hsa")#物种

可视化，其他的可视化调整可自行学习函数相关参数：

#结果可视化
data_KEGG <- pairwise_termsim(data_KEGG)#要做分组图，需要先运行这个函数
emapplot(data_KEGG, showCategory=8,legend_n=2)+
  scale_fill_manual(values = dittoColors())#修改填充颜色

image.png

关于富集分析的内容就补充到这里了，觉得分享有用的点个赞再走呗！