GO富集可视化：clusterProfiler GOplot

生物_医药_科研 2020-04-13

展开全文

现在的社会正处于大数据时代，生物学研究领域也迎来了大组学时代。无论是植物。动物，还是微生物，都有许多物种被测序。在后续的个性化分析过程中，通常会得到庞大的基因数据，如何选择一个高效的生物功能分析方法显得尤其重要。基因的功能富集，可以让你快速地了解所获得的基因具有怎样的生物功能，它们主要集中在哪些生物通路。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) #查看数据      V11 ADCY102  ADCY23  ADCY34  ADCY45  ADCY56  ADCY6gene=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转换为entrzidego <- 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   9ego2 <- 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'

###构建基因矩阵，并随机生成logFCgenedata = 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 itGOBubble(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'))