#GO柱状图绘制教程
#激活r包
library(GOplot)
  1. 数据准备与说明
    1.1 查看教程中给的数据
#查看教程中的数据
write.csv(EC$david,"EC-david.csv")
write.csv(EC$genelist,"EC-genelist.csv")

数据说明

R语言KEGG富集气泡图怎么按照P值升序 go富集柱状图r语言_数据集

EC$david为GO富集分析结果,包括174条显著富集的GOterm,第一列为序号,第二列为GO富集分析的类型,第三列为GOterm的ID,第四列为GOterm的名称,第五列使富集到这个term上的具体的基因(全大写),第六列使富集的显著性(校正P值)。

R语言KEGG富集气泡图怎么按照P值升序 go富集柱状图r语言_r语言_02

数据的EC$genelist为筛选出的差异表达的数据集,包括校正p值小于0.05的2040条基因,第一列为序号,第二列为基因名(首字母大写),第三列为差异倍数,第四列暂不清楚,第五列暂不清楚,第六列为P值,第七列为校正P值,第八列暂不清楚。

#生成绘图对象
circ <- circle_dat(EC$david, EC$genelist)
#查看一下生成的用于绘图的数据格式
write.csv(circ,"circ-david+genelist.csv")

R语言KEGG富集气泡图怎么按照P值升序 go富集柱状图r语言_r语言_03

可以看到,circle_dat函数将david和genelist两个数据集进行了整合,将david数据的每个富集到goterm里面的基因拆开了,然后生成了counts列,为每个基因更具genelist里的信息添加了差异倍数、校正p值以及自动计算出zscore值(具体计算方法在帮助文档里)。前面genelist里面那些不清楚的数据用不到。总共生成了9804行,这是因为有很多基因富集到了不同的term中。

#生成一个简单的barplot
GOBar(subset(circ, category == 'BP'))#只生成BP的富集结果

R语言KEGG富集气泡图怎么按照P值升序 go富集柱状图r语言_ci_04

##根据术语的类别来修饰barplot 
GOBar(circ, display = 'multiple')

R语言KEGG富集气泡图怎么按照P值升序 go富集柱状图r语言_ci_05

# 添加标题和改变颜色
GOBar(circ, display = 'multiple', title = 'Z-score coloured barplot', zsc.col = c('yellow', 'black', 'cyan'))

R语言KEGG富集气泡图怎么按照P值升序 go富集柱状图r语言_数据_06


更改完颜色后反而觉得变丑了,GOterm的名字也不清楚,具体的颜色可以更改代码个性化一下。好了,下面上一下自己的数据。

先按照上面数据的格式把表头改好(懒得指定了免得报错),然后不需要的数据给删掉(不删应该也行)。

然后运行上面代码,发现出来的图是这样的。

R语言KEGG富集气泡图怎么按照P值升序 go富集柱状图r语言_ci_07

很奇怪,很多富集到的通路直接消失了,因为这个柱状图是以校正p值给的y轴,z-score给的x轴,所以先关注数据集中这两个值有没有问题。

R语言KEGG富集气泡图怎么按照P值升序 go富集柱状图r语言_ci_08


发现一个term里面只要有一个基因的差异倍数的值缺失就会导致整个term的Zscore值缺失。查看一下zscore值得计算方法。

R语言KEGG富集气泡图怎么按照P值升序 go富集柱状图r语言_r语言_09


这个值是由上调基因个数减下调基因个数再除总数的开根。我们可以看到数据集中LBSP基因得差异倍数是缺失的,这可能是由于实验组或对照组该基因完全不表达使得计算差异倍数时报错。进而导致这一系列的报错。为了规避这种错误,可只计算有表达的那一组的表达量的对数。(或者直接给±1就好了,因为本身这个公式就只计算上调或者下调的个数)

还有一个问题是,david默认的是以P小于0.1进行富集的,也就是说,校正p值取对数直接约等于没有了。

问题解决,把校正p值小于0.05的term删了

R语言KEGG富集气泡图怎么按照P值升序 go富集柱状图r语言_r语言_10