这篇帖子其实是更新、补充、解决一下问题的。我们号写过很多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
关于富集分析的内容就补充到这里了，觉得分享有用的点个赞再走呗！
网友评论