rna_seqR语言可视化之美RNA-Seq

TCGA|GEO可视化分析第2篇---KEGG|GO分析大全

2020-04-19  本文已影响0人  沉迷工作的我

导读:之前讲了如何用筛选出的差异基因做做相关性分析。那今天我和各位小伙伴深入的讲一下:(1)如何用clusterProfiler做KEGG|GO富集条形图,气泡图;(2)如何用enrichplot做gene-GO terms|gene-KEGG pathways网络图;(3)如何用GOplot绘制gene-GO terms|gene-KEGG pathways和弦图;(4)如何用heatplot绘制gene-GO terms|gene-KEGG pathways关系瀑布图。


正文:GO,KEGG富集分析条形图(barplot)和气泡图(dotplot)都只显示最显著的富集项,而用户如果想知道哪些基因与这些显著项有关,该怎么办呢?cnetplot就是将基因GO terms或KEGG pathways之间的联系描述为一个网络,体现gene-GO terms|gene-KEGG pathways的关系趋势。瀑布图(heatplot)与cnetplot类似,它将gene-GO terms|gene-KEGG pathways关系显示为热图。因为gene-GO terms|gene-KEGG pathway网络图可能会过于复杂,热图可以简化结果,更容易识别表达模式。

Part1:KEGG|GO富集条形图,气泡图

library(openxlsx)

library(clusterProfiler)

library(ggplot2)

library(enrichplot)

library(GOplot)

library(DOSE)

library(stringr)

先把需要用到的包加载好,没有安装的依赖包一起安装一下。安装好所有的包之后,下面我们开始读取数据。

setwd("E:\\Bioinfo_analysis\\scripts\\corr\\corr_batch") #设置工作路径

fr<-read.xlsx('DEGs.xlsx',rowNames = F,colNames = T)#数据格式如下

gene<-bitr(fr$Genes,fromType = 'SYMBOL',toType = 'ENTREZID',OrgDb = 'org.Hs.eg.db') #基因名ID转换,把基因名转换成ENTREZID

GO<-enrichGO(

  gene$ENTREZID,

  OrgDb = 'org.Hs.eg.db',

  keyType = "ENTREZID",

  ont = "ALL",

  pvalueCutoff = 0.05,

  pAdjustMethod = "BH",

  qvalueCutoff = 0.05,

  minGSSize = 10,

  maxGSSize = 500,

  readable = TRUE

)#GO富集

KEGG<-enrichKEGG(

  gene$ENTREZID,

  organism = "hsa",#我用到是数据是人的组织数据,所以这里选择‘hsa’

  keyType = "kegg",

  pvalueCutoff = 0.05,

  pAdjustMethod = "BH",

  universe,

  minGSSize = 5,

  maxGSSize = 500,

  qvalueCutoff = 0.2,

  use_internal_data = FALSE

)#KEGG富集

barplot(GO, split="ONTOLOGY")+facet_grid(ONTOLOGY~., scale="free") #GO聚类条形图,fig1

fig1.BP:生物过程,CC:细胞组成,MF:分子功能。

dotplot(GO, split="ONTOLOGY")+facet_grid(ONTOLOGY~., scale="free") #GO聚类气泡图,fig2

Fig2

barplot(KEGG,showCategory = 40,title = 'KEGG Pathway') #KEGG聚类条形图,fig3

fig3

dotplot(KEGG) #KEGG聚类气泡图,fig4

fig4

Part2:gene-GO terms|gene-KEGG pathways网络图

其实以上都是KEGG和GO相关的常规分析,GO,KEGG富集分析条形图(barplot)和气泡图(dotplot)都只显示最显著的富集项,而用户如果想知道哪些基因与这些显著项有关,该怎么办呢?下面我们再升级一下。

enrichplot::cnetplot(GO,circular=TRUE,colorEdge = TRUE)#GO通路-基因网络图,fig5

fig5

enrichplot::cnetplot(KEGG,circular=TRUE,colorEdge = TRUE)#KEGG通路-基因网络图,fig6

fig6

Fig5和fig6体现的就是KEGG和GO每个富集项中都有哪些基因,我们由此也可以推测一些基因的功能。

Part3:gene-GO terms|gene-KEGG pathways关系瀑布图

我们筛选的差异基因还是比较少的,这样的网路图看起来还是很容易区分不同的表达模式的,不过当差异基因过多的时候,这样的网络图就会显得比较乱,可能过于复杂,那我们应该怎么办呢?不用慌,我们可以用瀑布图来展示。瀑布图(heatplot)与cnetplot类似,它将gene-GO terms|gene-KEGG pathways关系显示为热图,热图可以简化结果,更容易识别表达模式。

enrichplot::heatplot(GO,showCategory = 50) #GO富集瀑布图,fig7

fig7

enrichplot::heatplot(KEGG,showCategory = 50) #kegg富集瀑布图,fig8

fig8

从Fig7和fig8中就可以更容易识别不同基因的表达模式。

Part4:gene-GO terms|gene-KEGG pathways和弦图

有些小伙伴说想让图能体现功能的同时更炫酷一点。OK,我们还可以再升级一下,此时用和弦图展示你的结果,定会让你的图片逼格升一个档次。好了,直接上代码。

GOplotIn<-GO[1:10,c(2,3,7,9)]#我们先提取GO富集结果的前10行,和提取ID,Description,p.adjust,GeneID四列。

GOplotIn$geneID <-str_replace_all(GOplotIn$geneID,'/',',')#把GeneID列中的’/’替换成‘,’

names(GOplotIn)<-c('ID','Term','adj_pval','Genes')#修改列名,后面和弦图绘制的时候需要这样的格式,不然会报错

GOplotIn$Category = "BP"#因为我们提取的前10列为BP,所以再加一列分类信息

#构建基因表达矩阵

genedata<-data.frame(ID=fr$Genes,logFC=fr$logFC)

circ<-GOplot::circle_dat(GOplotIn,genedata) #GOplot导入数据格式整理

chord<-chord_dat(data = circ,genes = genedata) #生成含有选定基因的数据框

好了,到此数据都整理好了,开始画和弦图。

GOChord( #GO富集和弦图,fig9

  data = chord,

  title = 'GOchord plot',

space = 0,#GO Term间距

  limit = c(1,1),

  gene.order = 'logFC',

  gene.space = 0.25,

  gene.size = 5,

lfc.col = c('red','white','blue'), #上下调基因颜色

  #ribbon.col = brewer.pal(length(GOplotIn$Term)),#GO Term colors

process.label = 10 #GO Term字体大小

)

Fig9是不是感觉很酷炫,撒花~~~~~~~~~~~~~~~~~~

还可以和弦图+表格一起展示GO富集结果,往下看:

GOCircle(circ) #和弦图+表格,fig10 

fig10

也可以聚类热图的形式展示,GO富集结果:

GOCluster(circ,GOplotIn$Term) #GO富集聚类图,fig11

fig11

以上就是GO和KEGG富集相关的分析,学会这些足以了,以后再也不用担心怎么去做GO,KEGG分析了。喜欢的小伙伴可以给个赞哦,有疑问和建议可以后方留言哦!

上一篇下一篇

猜你喜欢

热点阅读