❝详情请联系作者：
❞ 这篇帖子其实是更新、补充、解决一下问题的。我们号写过很多GO、KEGG富集分析的可视化，数不胜数，可以在公众号检索“富集”了解更多。我们演示的时候都是直接提供了富集的结果文件，一般演示为了图方便，也是利用在线工具cytoscape做的。结果一伙伴最近提问有做过GO、KEGG富集的R语言帖子么，突然发现这样的内容还没有正经写过，所以这里补充一下。 首先我们做一下单独的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./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/aPEARlibrary(aPEAR)#相关阈值设定清查阅读 ?enrichmentNetwork#enrichmentNetwork返回的是ggplot对象，所以可以用ggplot2主题修饰enrichmentNetwork(sig_gene_KEGG[[1]], colorBy = 'pvalue', colorType = 'pval')+ scale_color_gradientn(colours = c("#B83D3D",'white','#1A5592'), name = "pvalue")
到这里我们这篇帖子还没结束，有两个原因。第一我们之前总是说一个函数，多组GO分析（例如：Clusterprofile多分组富集分析及可视化），但是没有说过KEGG，有小伙伴去做的时候出错，其实参数里面选择KEGG，设置对就可以进行这个分析。第二个问题是有小伙伴发来图让复现，是富集结果的展示，乍一看很复杂，既是网络图，又是多组的，其实很简单，clusterProfiler多组富集分析和enrichplot早就解决了这个问题。所以这两个问题我们归为一个进行解决。 
（reference：A single-cell transcriptomic atlas of exercise-induced antiinflammatory and geroprotective effects across the body） 我们还是利用之前的分组基因的文件，进行多组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转化为IDGene_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())#修改填充颜色
关于富集分析的内容就补充到这里了，觉得分享有用的点个赞再走呗！
|
