GEO生信分析流程小教程收藏

富集分析后的结果如何画chord弦状图

2020-06-17  本文已影响0人  小梦游仙境

Chord Diagram弦状图-用来画富集分析的结果

前面看的文献中,有一个很漂亮的圈图,https://mp.weixin.qq.com/s/NIqSTvJ916ZGk_dwsRlH5w,就是kegg和go分析后的基因和通路之间的关系图,用cnetplot函数画出来的图也是一样的含义,不过这个圈图看起来就真的很好看,就想画,主要其实颜色是可以选择的,用https://www.sioe.cn/yingyong/yanse-rgb-16/中的颜色就可以选择,所以可以画出来更好看的图。不过后面在数据变换的时候,就还是没转过来可以直接用最后面画图用的数据,就还在思考如何一步步构建包中的中间数据从而得到足后的数据,后来老大告诉我直接构建后面的数据就可以了。

弦状图这种类型的图表可视化了实体之间的相互关系。实体之间的连接用于显示它们共享某些相同的东西。节点沿圆排列,通过使用圆弧将点之间的关系连接在一起。值被分配给每个连接,它由每个弧的大小成比例地表示。颜色可以用来把数据分组到不同的类别,这有助于进行比较和区分组。不过当有太多的连接显示时,可能就会出现混乱。

这使得弦图非常适合于比较数据集内或不同数据组之间的相似性。

使用R语言中的绘图包,比如 circlize, RCircos, CIRCUS, and OmicCircos等来绘制。

一张显示人口流动的图如下,来自http://download.gsb.bund.de/BIB/global_flow/

image-20200617112544546

GOplot包可用于生物数据的可视化。但是要注意该包不能用于执行这些分析,只能把分析结果进行可视化

这个包的用法参考:http://wencke.github.io/。主要是可以把富集分析后的结果进行圈图展示。链接中的讲解代码很长,提炼后的代码如下

david <- EC$david
genelist <- EC$genelist
circ <- circle_dat(EC$david, EC$genelist)
chord <- chord_dat(data = circ, genes = EC$genes, process = EC$process)
GOChord(chord, space = 0.02, gene.order = 'logFC', gene.space = 0.25, gene.size = 5)

画出的图如下

image-20200617120210980

EC是这个GOplot包中的实例数据,为了能用这个包中的作图函数,就先了解这个包中的数据。那么david和genelist的数据长下面这样。

image-20200617114924282

生成的中间数据chord长下面这样。

image-20200617115516739

最后用GOChord函数画图的数据chord的数据格式如下图所示。这是我们要把自己的数据变成这样的数据的样子。就是下面的chord数据格式。变成下面这样就可以画圈图了。

image-20200617114629848

上面最后画图的输入数据函数就是chord,最后的logFC是从前面的差异分析的结果得来的,也就是genelist中的logFC。然后chord数据除了logFC以外,前面的数据就是如果某个基因在通路中,就是1,否则就是0。这中数据就是最前面代码中的chord_dat函数实现的,但是问题是这个chord 数据是从circ得到的,而circ是从david和genelist得到的,david数据是富集分析的结果,genelist是差异分析的结果,现在我没有差异分析的结果,我只是有前面比如通过cox回归筛选得到的基因名,怎么办呢?

问题解决的办法就是,直接制作数据格式到chord这样的数据,把logFC变为0就可以了。

load(file = 'sur_logKM_cox.Rdata')
cox=rownames(cox_results[cox_results[,4]<0.01,])
log=names(log_rank_p[log_rank_p<0.01])
surGenes=intersect(cox,log)

上面得到的surGenes如下:就是一堆字符串啦

image-20200617121812423

但是用老大出品的这个AnnoProbe进行ID转换,得到了如下图所示非常方便后续进行操作的表格


library(AnnoProbe) #用老大出品的这个AnnoProbe进行ID转换,
surGenes=annoGene(surGenes,ID_type = 'SYMBOL','human') 


tail(sort(table(surGenes$biotypes)))
head(surGenes)

image-20200617122118975

接下来富集分析

library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
df <- bitr(unique(surGenes$SYMBOL), fromType = "SYMBOL",
           toType = c( "ENTREZID"),
           OrgDb = org.Hs.eg.db)#富集分析要求数据数据的是ENTREZID,要从前面得到的SYMBOL用bitr函数得到ENTREZID
head(df)   

gene_up=df$ENTREZID

#下面是为了画图做的插曲
#将bp的结果保存
bp_go <- enrichGO(gene_up, 
                  OrgDb = "org.Hs.eg.db", 
                  ont="BP",
                  pvalueCutoff  = 0.001,
                  qvalueCutoff  = 0.001) 

enrichgo=DOSE::setReadable(bp_go, OrgDb='org.Hs.eg.db',keyType='ENTREZID')#setReadable函数直接进行ID转换了,从ENTREZID再得到symbol

mm <- enrichgo@result[1:5,c(2,7,8)]

得到的mm数据如下(mm是我名字的缩写,哈哈)

image-20200617122407370
tmp=do.call(rbind,
        apply(mm, 1,function(x){
          data.frame(go=x[1],
                     gene=strsplit(x[3],'/')[[1]])
        })
)

得到的tmp数据格式如下

image-20200617123212419

接下来就是将这个长形数据用reshape包中的dcast函数变为宽形数据,如下

library(reshape2)
tmp2=dcast(tmp,go~gene)
tmp2[is.na(tmp2)]=0
rownames(tmp2)=tmp2[,1]
tmp2=tmp2[,-1]
tmp2=t(tmp2)
tmp2[tmp2!=0]=1
tmp2=as.data.frame(tmp2)
tmp2$logFC=0
cg=rownames(tmp2)
tmp2=apply(tmp2,2,as.numeric)
rownames(tmp2)=cg

得到的tmp2数据如下,这就是已经做好了和前面的实例数据chord长得一摸一样的数据了。

image-20200617123041839

接下来画图

GOChord(tmp2, space = 0.02, gene.order = 'logFC', gene.space = 0.25, gene.size = 5)
image-20200617123440677

区别是我们没有用变化倍数进行聚类,所以没有logFC的渐变的变化,不过已经超级好看了!

如果想有倍数变化的结果,那么就是一开始就是有差异分析的结果,就可以下面的代码,得到的就是本文第二张图片了。

# Create the plot
chord <- chord_dat(data = circ, genes = EC$genes, process = EC$process)
GOChord(chord, space = 0.02, gene.order = 'logFC', gene.space = 0.25, gene.size = 5)

参考的代码全套在http://wencke.github.io/,也可以搜索到翻译过来的中文版的教程。

上一篇下一篇

猜你喜欢

热点阅读