[组学分析系列的教程]，前期还是整理一下前期的教程。前期自己分享了精美图形的绘制。但是，距离现在很久了，不知道有多少童鞋还会使用或回去复习。从今天开始，前期我们一起学习这个教程，我们的目标是：至少可以满足现在的自己使用！！
?
本期分享的是：功能富集分析。
?
---
?
我们童鞋最关心的一件事：如何获得绘图代码和数据！！
途径一：查看我们以前的教程
途径二：给本文打赏，小杜看到后立即将本文的实例数据和代码发给你！
?
---
?
从分享自己的学习教程至今快有小一年了，我的账号都是没有使用超能力推广的，因此关注的人还是处于“小白”阶段，与我自己对自己的定义是一样的。自己的更新是不定时的，也许你会持续看到我更新，也许你很久看不到我更新文章，这都是正常状态，自己一个人的精力有限，请谅解。也希望大家可以投稿，分享自己的学习笔记。
此外，有童鞋问：是否有“粉丝群”，我的天啊！我也会有一天可以建立自己的“群”！但是，现在我暂时是不会建群的！因为，自己能力有限。不能很好的帮助大家，以及建立这些群管理起来也是比较费时间的。
但是，但是，如果你愿意，可以加小杜的微信，【前提：需转账50元（不介意更多）】。目前，只要是小杜的微信好友，都可以直接向小杜索要每期的教程的实例数据和代码（这部分也是后期小杜“粉丝群”的首选人员，哈哈哈哈........）[微信可在公众号界面获得]。
温馨提示：这个请根据自己的意愿选择，不要盲目的添加。小杜也是生信小白，也许无法提供给你想要的帮助，请理智消费。
?
---
GO功能富集热图
GO网络图
网络图
通路图
---
图形绘
## clusterProfiler功能富集
## author:小杜的生信筆記
setwd("D:\\小杜的生信筆記\\clusterProfiler功能富集")
##导入相关的R包
library(stringr)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(topGO)
library(enrichplot)
## 导入数据
df <- read.csv("DEGs_ID.csv", header = F)
head(df)
dim(df)
## 查看数据类型
class(df)
## 转换数据类型
df02 <- as.character(df$V1)
head(df02)
class(df02)
## symbol To entrez ID
DEG.gene_symbol = df02
DEG.entrez_id = mapIds(x = org.Hs.eg.db,
keys = DEG.gene_symbol,
keytype = "SYMBOL",
column = "ENTREZID")
gene = bitr(DEG.gene_symbol, fromType="SYMBOL",
toType="ENTREZID",
OrgDb="org.Hs.eg.db")
#write.csv(gene, "SYMBOL_WNTREZID.list.csv")
head(gene)
## GO 富集分析
ego <- enrichGO(
gene = gene$ENTREZID, ## ENTREZID
keyType = "ENTREZID", ## 输入基因类型
OrgDb = org.Hs.eg.db, ## 导入背景基因
ont = "all", ## GO的种类，BP，CC， MF
pAdjustMethod = "BH", ## 矫正的类型
pvalueCutoff = 0.05, ## P值过滤值
readable = TRUE)
head(ego)
### GO 注释柱状图
pdf(file = "GO_bar.pdf", width = 10, height = 9)
barplot(ego,
drop = TRUE,
showCategory = 10, ## 显示前10个GO term
split="ONTOLOGY") +
scale_y_discrete(labels=function(x) str_wrap(x, width=80))+
facet_grid(ONTOLOGY~., scale='free')
dev.off()
## GO 气泡图
pdf(file = "GO_dotplot.pdf", width = 8, height = 8)
dotplot(ego, showCategory = 10)
dev.off()
# 导出GO数据
write.csv(ego, "GO.output.csv")
##------------------------------------------------------------------------------
## GO网络图
BP_ego <- enrichGO(
gene = gene$ENTREZID,
keyType = "ENTREZID",
OrgDb = org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
readable = TRUE)
## BP 富集网络图
pdf("BP_网络图.pdf", width = 10, height = 10)
plotGOgraph(BP_ego)
dev.off()
### 使用simplify()函数过滤
BP_ego.filt <- simplify(BP_ego,
cutoff = 0.7,
by = "p.adjust",
select_fun = min,
measure = "Wang",
semData = NULL)
head(BP_ego.filt)
pdf("BP_网络图.filt.pdf", width = 10, height = 10)
plotGOgraph(BP_ego.filt)
dev.off()
## 好看的网络图
BP_ego <- enrichGO(
gene = gene$ENTREZID,
keyType = "ENTREZID",
OrgDb = org.Hs.eg.db,
ont = "BP",
readable = TRUE)
## visualization
BP_ego <- pairwise_termsim(BP_ego)
pdf("好看的GO网络图.test.pdf", width = 9, height = 10)
emapplot(BP_ego, cex_label_category=.8, cex_line=.5) +
scale_fill_continuous(low = "#e06663", high = "#327eba", name = "p.adjust",
guide = guide_colorbar(reverse = TRUE, order=2.5), trans='log10')
print(p4)
dev.off()
## use simplify to remove redundant terms
BP_ego2 <- simplify(BP_ego, cutoff=0.8, by="p.adjust", select_fun=min)
BP_ego2 <- pairwise_termsim(BP_ego2)
pdf("DK.filt04.pdf", width = 8, height = 8)
emapplot(BP_ego2, cex_label_category=.8, cex_line=.5) +
scale_fill_continuous(low = "#8da0cb", high = "#1b9e77", name = "p.adjust",
guide = guide_colorbar(reverse = TRUE, order=1.4), trans='log10')
dev.off()
## -----------------------------------------------------------------------------
## KEGG 富集分析
ekegg <- enrichKEGG(
gene = gene$ENTREZID,
keyType = "kegg",
organism = 'hsa',
pvalueCutoff = 0.05,
pAdjustMethod = "BH")
## KEGG 柱状图
pdf(file = "KEGG_bar.pdf", width = 9, height = 8)
barplot(ekegg, showCategory = 10)
dev.off()
## KEGG 气泡图
pdf(file = "KEGG_dotplot.pdf", width = 9, height = 8)
dotplot(ekegg, showCategory = 10)
dev.off()
## 导出KEGG富集分析数据
write.csv(ekegg, "KEGG.output.csv")
### -----------------
## 气泡图 更改颜色
dotplot(ekegg, showCategory = 10)+
scale_color_gradient(low ="#e06663", high = "#327eba", name = "p.adjust" )
##
library(ggplot2)#柱状图和点状图
library(stringr)#基因ID转换
library(enrichplot)#GO,KEGG,GSEA
library(clusterProfiler)#GO,KEGG,GSEA
library(GOplot)#弦图，弦表图，系统聚类图
library(DOSE)
library(ggnewscale)
library(topGO)#绘制通路网络图
library(circlize)#绘制富集分析圈图
library(ComplexHeatmap)#绘制图例
###
GO <- enrichGO(
gene = gene$ENTREZID, ## ENTREZID
keyType = "ENTREZID", ## 输入基因类型
OrgDb = org.Hs.eg.db, ## 导入背景基因
ont = "ALL", ## GO的种类，BP，CC， MF
pAdjustMethod = "BH", ## 矫正的类型
pvalueCutoff = 0.05, ## P值过滤值
qvalueCutoff = 0.05,#设定q值阈值
readable = TRUE)
pdf("GO网络图.pdf", width = 20, height = 10)
enrichplot::cnetplot(GO,circular=FALSE,colorEdge = TRUE)#基因-通路关联网络图
dev.off()
pdf("KEGG网络图.pdf", width = 10, height = 10)
enrichplot::cnetplot(ekegg,circular=FALSE,colorEdge = TRUE)#circluar为指定是否环化，基因过多时建议设置为FALSE
dev.off()
## 也可以以热图形式展现关联关系:
pdf("GO富集热图.pdf", width = 20, height = 15)
enrichplot::heatplot(GO,showCategory = 20)#基因-通路关联热图
dev.off()
pdf("KEGG富集热图.pdf", width = 20, height = 15)
enrichplot::heatplot(ekegg,showCategory = 20)
dev.off()
## 富集到的功能集/通路集之间的关联网络图：
GO2 <- pairwise_termsim(GO)
KEGG2 <- pairwise_termsim(ekegg)
enrichplot::emapplot(GO2,showCategory = 50, color = "p.adjust", layout = "kk")#通路间关联网络图
enrichplot::emapplot(KEGG2,showCategory =50, color = "p.adjust", layout = "kk")
## kegg 选择通路进行展示:
browseKEGG(ekegg,"hsa04640")#选择其中的hsa05166通路进行展示
## GO富集弦图:
write.table(GO$ONTOLOGY, "GO_ONTOLOGYs.txt", #将所有GO富集到的基因集所对应的类型写入本地文件从而得到BP/CC/MF各自的起始位置如我的数据里是1，2103，2410
append = FALSE, quote = TRUE, sep = " ",
eol = "\n", na = "NA", dec = ".", row.names = TRUE,
col.names = TRUE, qmethod = c("escape", "double"),
fileEncoding = "")
GOplotIn_BP<-GO[1:10,c(2,3,7,9)] #提取GO富集BP的前10行,提取ID,Description,p.adjust,GeneID四列
GOplotIn_CC<-GO[642:710,c(2,3,7,9)]#提取GO富集CC的前10行,提取ID,Description,p.adjust,GeneID四列
GOplotIn_MF<-GO[711:775,c(2,3,7,9)]#提取GO富集MF的前10行,提取ID,Description,p.adjust,GeneID四列
GOplotIn_BP$geneID <-str_replace_all(GOplotIn_BP$geneID,'/',',') #把GeneID列中的’/’替换成‘,’
GOplotIn_CC$geneID <-str_replace_all(GOplotIn_CC$geneID,'/',',')
GOplotIn_MF$geneID <-str_replace_all(GOplotIn_MF$geneID,'/',',')
names(GOplotIn_BP)<-c('ID','Term','adj_pval','Genes')#修改列名,后面弦图绘制的时候需要这样的格式
names(GOplotIn_CC)<-c('ID','Term','adj_pval','Genes')
names(GOplotIn_MF)<-c('ID','Term','adj_pval','Genes')
GOplotIn_BP$Category = "BP"#分类信息
GOplotIn_CC$Category = "CC"
GOplotIn_MF$Category = "MF"
circ_BP<-GOplot::circle_dat(GOplotIn_BP,genedata) #GOplot导入数据格式整理
circ_CC<-GOplot::circle_dat(GOplotIn_CC,genedata)
circ_MF<-GOplot::circle_dat(GOplotIn_MF,genedata)
chord_BP<-chord_dat(data = circ_BP,genes = genedata) #生成含有选定基因的数据框
chord_CC<-chord_dat(data = circ_CC,genes = genedata)
chord_MF<-chord_dat(data = circ_MF,genes = genedata)
GOChord(data = chord_BP,#弦图
title = 'GO-Biological Process',space = 0.01,#GO Term间距
limit = c(1,1),gene.order = 'logFC',gene.space = 0.25,gene.size = 5,
lfc.col = c('red','white','blue'), #上下调基因颜色
process.label = 10) #GO Term字体大小
GOChord(data = chord_CC,title = 'GO-Cellular Component',space = 0.01,
limit = c(1,1),gene.order = 'logFC',gene.space = 0.25,gene.size = 5,
lfc.col = c('red','white','blue'),
process.label = 10)
GOChord(data = chord_MF,title = 'GO-Mollecular Function',space = 0.01,
limit = c(1,1),gene.order = 'logFC',gene.space = 0.25,gene.size = 5,
lfc.col = c('red','white','blue'),
process.label = 10)
##
chord<-chord_dat(data = circ_BP,genes = df) #生成含有选定基因的数据框
GOCluster(circ_BP,GOplotIn_BP$Term) #系统聚类图---
往期文章
---
小杜的生信筆記 ，主要发表或收录生物信息学的教程，以及基于R的分析和可视化（包括数据分析，图形绘制等）；分享感兴趣的文献和学习资料！
