现在的社会正处于大数据时代，生物学研究领域也迎来了大组学时代。无论是植物。动物，还是微生物，都有许多物种被测序。在后续的个性化分析过程中，通常会得到庞大的基因数据，如何选择一个高效的生物功能分析方法显得尤其重要。基因的功能富集，可以让你快速地了解所获得的基因具有怎样的生物功能，它们主要集中在哪些生物通路。clusterProfiler、GOplot这两个工具可以说是在众多GO富集分析软件中脱颖而出，它们不仅能够支持多个数据库（GO/KEGG/DOSE/MSigDb/Reactome），还能绘出高大上的CNS级别富集图。那么，今天主要就是介绍如何结合这两个软件，让你的一组gene symbol秒变高富帅。

01. 准备数据

经过上游的生信分析我们会获得许多具有生物学意义的gene set，可以是差异表达基因，也可是正选择基因或者加速进化基因。通常，只要具有这些基因的gene symbol或者是geneid，都可以利用该软件进行分析。

本次所使用的数据主要是gene symbol，如下图所示

02. 读入数据

setwd("C:/Users/lenovo/Desktop/RERconverge-master/R")  # 设置工作路径
geneset <- read.table(file = "genelist",header = F) #读取数据
> head(geneset) #查看数据
      V1
1 ADCY10
2  ADCY2
3  ADCY3
4  ADCY4
5  ADCY5
6  ADCY6
gene=as.character(geneset[,1]) #转换字符格式
> head(gene)
[1] "ADCY10" "ADCY2"  "ADCY3"  "ADCY4" 
[5] "ADCY5"  "ADCY6"

03. 使用clusterProfiler进行富集分析

library(org.Hs.eg.db)
library(clusterProfiler)
library(GOplot)
library(ggplot2)
library(stringr)

gene <- bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db") # 将symbol转换为entrzid

ego <- enrichGO(gene=gene$ENTREZID,
                OrgDb='org.Hs.eg.db',
                ont= "BP",
                pAdjustMethod="BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
                 readable = TRUE)
> dim(ego)
[1] 402   9
ego2 <- simplify(ego,cutoff=0.7,by="p.adjust",select_fun=min)  #去除冗余，可以调整cutoff值
> dim(ego2)
[1] 144   9

heatlpot(ego2)

cnetplot(ego2,categorySize="pvalue",circular = TRUE, colorEdge = TRUE)

# KEGG富集
kk <- enrichKEGG(gene$ENTREZID, organism="hsa",keyType = "kegg",pvalueCutoff=0.01,pAdjustMethod="BH",qvalueCutoff=0.05)
barplot(kk, showCategory=20,title="Enrichment KEGG")

04. 整理GOplot输入文件格式

###ego3 <- data.frame(ego2) ### 转换为数据矩阵
GO <- ego[1:9,c(1,2,8,6)]  ### 选取富集结果前9个，并提取GO id，description， adj.p以及基因ID列

GO$geneID <- str_replace_all(GO$geneID,"/",",") ### 修改geneID这一列的分隔符号
names(GO)=c("ID","Term","Genes","adj_pval")
GO$Category = "BP"

###构建基因矩阵，并随机生成logFC
genedata = data.frame(ID=geneset$V1,logFC=rnorm(length(geneset$V1),mean = 0,sd=2))
circ <- circle_dat(GO,genedata)

05. 使用GOplot绘图

# 01. 和弦图
chord <- chord_dat(data=circ, genes = genedata)  # 生成带有选定基因列表的矩阵
chord <- chord_dat(data=circ, process = GO$Term) #生成带有选定GO term的列表矩阵

chord <- chord_dat(data=circ, genes=genedata,process = GO$Term) # 构建数据
GOChord (         data=chord,
                  title = "GOChord plot", 
                  space = 0.02, # go term处间隔大小
                  limit = c(3,5) #第一个值是至少分配给一个基因的go term数目，第二个数值是至少分配给一个Go term的基因数
                  gene.order ='logFC',gene.space=0.25,gene.size=10,
                  lfc.col = c('firebrick3','white','royalblue3') #上下调基因颜色
                  ribbon.col=brewer.pal(length(GO$Term)),"set3" #GO term颜色
                  process.label = 18 # Go terms字体大小
)


# 02. 条形图
GOBar(circ,display = 'multiple')

# 03. 气泡图
# 要添加标题，请更改圆圈的颜色，对图进行构图，并更改标签阈值，请使用以下参数：
GOBubble(circ, title = 'Bubble plot', colour = c('orange', 'darkred', 'gold'), display = 'multiple', labels = 3)
# 对于构面图，还可以通过将bg.col设置为TRUE，根据显示的类别为面板的背景着色:
GOBubble(circ, title = 'Bubble plot with background colour', display = 'multiple', bg.col = T, labels = 3)
# 软件包的更新版本中包含一个新函数reduce_overlap ，以减少冗余项的数量. 到目前为止，已实现的方法非常简单+缓慢，需要进一步完善. 但是，通过减少冗余项的数量，可以像气泡图一样显着提高图的可读性. 该函数删除基因重叠大于或等于设定阈值的所有术语. 该函数在不考虑GO层次结构的情况下，每组代表一个术语:
# Reduce redundant terms with a gene overlap >= 0.75...
reduced_circ <- reduce_overlap(circ, overlap = 0.75)
# ...and plot it
GOBubble(reduced_circ, labels = 2.8)

# 04. 圈图
GOCircle(circ)

# 外圈显示了分配基因的logFC的每个项的散点图. 默认情况下，红色圆圈显示上调，蓝色圆圈显示下调. 可以使用参数lfc.col更改颜色. 因此，更容易理解，为什么在某些情况下，高度有效的术语的z得分接近于零. Z分数为零并不意味着该术语不重要. 至少没有，只要其显着丰富即可. 它只是表明z分数是一个粗略的度量，因为显然分数并未考虑过程中单个基因的功能水平和激活依赖性. 您可以使用各种参数来更改图的布局，请参阅？ GOCirlce.nsub参数需要更多说明才能明智地使用. 首先，它可以是数字或字符向量. 如果它是一个字符向量，则它包含要显示的进程的ID或术语说明（未显示输出）。
IDs <- c('GO:0007507', 'GO:0001568', 'GO:0001944', 'GO:0048729', 'GO:0048514', 'GO:0005886', 'GO:0008092', 'GO:0008047')
GOCircle(circ, nsub = IDs)

# 05. 热图
GOHeat(chord, nlfc = 1, fill.col = c('red', 'yellow', 'green'))

# 06. 聚类图
GOCluster(circ, EC$process, clust.by = 'logFC', term.width = 2)

OCluster(circ, EC$process, clust.by = 'term', lfc.col = c('darkgoldenrod1', 'black', 'cyan1'))

# 07. 韦恩图

l1 <- subset(circ, term == 'heart development', c(genes,logFC))
l2 <- subset(circ, term == 'plasma membrane', c(genes,logFC))
l3 <- subset(circ, term == 'tissue morphogenesis', c(genes,logFC))
GOVenn(l1,l2,l3, label = c('heart development', 'plasma membrane', 'tissue morphogenesis'))

参考链接：

• https://yulab-smu.github.io/clusterProfiler-book/index.html
• http://wencke.github.io/
• http://www.360doc.com/content/19/0416/12/52645714_829160785.shtml
• https://www.jianshu.com/p/77246ff36214
• https://bioinfohome.com/index.php/2019/07/08/goplot-01/
• http://www.sohu.com/a/339298858_120098972
• http://wap.sciencenet.cn/blog-118204-1192033.html?mobile=1
• http://cran.r-project.org.icopy.site/web/packages/GOplot/vignettes/GOplot_vignette.html
• https://www.jianshu.com/p/16fd75aff9d4?utm_campaign=hugo&utm_medium=reader_share&utm_content=note&utm_source=qq