专栏名称: 生信人
共同学习生物信息学知识,共同探究生物奥秘。
目录
相关文章推荐
华大集团BGI  ·  援外培训 | ... ·  2 天前  
华大集团BGI  ·  2024国际衰老与再生研讨会在华大时空中心举行 ·  3 天前  
生物探索  ·  Nature | ... ·  5 天前  
生物制品圈  ·  生物制品圈招新啦,期待您的加盟! ·  6 天前  
BioArt  ·  Cell Metab | ... ·  6 天前  
51好读  ›  专栏  ›  生信人

用小提琴图来可视化基因差异表达之美化过程

生信人  · 公众号  · 生物  · 2017-07-06 07:24

正文

为了表示基因在样本中的差异,对于许多个基因的话用火山图感觉挺高大上的,那么对于少数几个基因的话一般用什么来展示其表达差异呢,有很多种办法,比如箱线图,柱状图等等,今天主要用小提琴图来可视化几个基因的差异表达。

首先安装R包

  1. source("https://bioconductor.org/biocLite.R")

  2. biocLite('vioplot')


然后我的数据是一个矩阵,行为样本列为基因如:

已经知道癌症样本和正常样本的索引

那么绘制小提琴图的简单办法如下:

  1. vioplot(sigData[tumor,1]

  2.    ,sigData[normal,1]

  3.    ,sigData[tumor,2]

  4.    ,sigData[normal,2]

  5.    ,sigData[tumor,3]

  6.    ,sigData[normal,3]

  7.    ,sigData[tumor,4]

  8.    ,sigData[normal,4]

  9.    ,names=rep(c("Tumor","Normal"),4),col="red")


从图中可以看出这四个基因其实在癌与癌旁是有差异的,但是好丑,很显然横轴要标上基因,每一个基因的癌与癌旁近一点,同时标记上p值最好了

那么怎么调整呢

首先需要调整的小提琴图的间距,让其看起来像四组,而不是八组,我们机智的设计了x轴,如代码所示:

  1. x=c(1:5)

  2. y=c(1:5)

  3. plot(x, y, xlim=c(0,10),xlab='',ylab='Gene Expression'

  4. ,main=''

  5. , ylim=c(min(sigData),max(sigData)+1),pch=21,col='white',xaxt="n")

  6. for(i in 1:4){

  7.  basal=sigData[tumor,i]

  8.  noBasal=sigData[normal,i]

  9.  vioplot(add = T,basal,at=3*(i-1),lty=1)

  10.  vioplot(add = T,noBasal,at=3*(i-1)+1,lty=1)

  11. }


从代码中可以看出我们将x轴设计成了0-10,然后再0,1,2,3,5,6,8,9处画上小提琴,这样就可以看出四组数据了,顺便将y轴的label加上

但是这个还是很丑,我们想将tumor和normal用颜色区分一下下,怎么办,如代码:

  1. x=c(1:5)

  2. y=c(1:5)

  3. plot(x, y, xlim=c(0,10),xlab='',ylab='Gene Expression'

  4. ,main=''

  5. , ylim=c(min(sigData),max(sigData)+1),pch=21,col='white',xaxt="n")

  6. for(i in 1:4){

  7.  basal=sigData[tumor,i]

  8.  noBasal=sigData[normal,i]

  9.  vioplot(add = T,col = 'blue',basal,at=3*(i-1),lty=1)

  10.  vioplot(add = T,col = 'red',noBasal,at=3*(i-1)+1,lty=1)

  11. }


如咱们所愿,颜色变了,但是p值木有加上,还是挺遗憾的,咱们使用‘Mann-Whitney’ test来检验差异显著性,然后把p值在图中,做法如下:

  1. x=c(1:5)

  2. y=c(1:5)

  3. plot(x, y, xlim=c(0,10),xlab='',ylab='Gene Expression'

  4. ,main=''

  5. , ylim=c(min(sigData),max(sigData)+1),pch=21,col='white',xaxt="n")

  6. for(i in 1:4){

  7.  basal=sigData[tumor,i]

  8.  noBasal=sigData[normal,i]

  9.  vioplot(add = T,col = 'blue',basal,at=3*(i-1),lty=1)

  10.  vioplot(add = T,col = 'red',noBasal,at=3*(i-1)+1,lty=1)

  11.  p=round(wilcox.test(basal,noBasal)$p.value,3)

  12.  mx=max(c(basal,noBasal))

  13.  lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx))

  14.  text(x=3*(i-    1)+0.5,y=mx+0.5,labels=ifelse(p==0,paste0("p<0.001"),paste0("p,p)),cex = 0.5)

  15. }



从图中可以看出p值已经加上去了,对比之前的代码主要多了四行代码

  1. p=round(wilcox.test(basal,noBasal)$p.value,3)#用来检验基因的差异显著性p值,并去3位有效数字

  2. mx=max(c(basal,noBasal))#找到新画上去的两个小提琴的最高的位置,为画上p值的位置坐准备

  3. lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx))#在画p值之前先画一条横跨在癌与癌旁的两个小提琴之间的横线,注意那个0.8,两个小提琴直接的间距是1,那么0.8表达离第二个小提琴的位置还差0.2

  4. text(x=3*(i-1)+0.5,y=mx+0.5,labels=ifelse(p==0,paste0("p<0.001") ,paste0("p,p)),cex = 0.5)#画上p值的文本啦,xy都是图中位置,仔细体会,label就是文本,cex设成0.5表示原来的50%

p值也画上去了,但是x轴一直没有,最后咱们来绘制x轴,首先分析一下八个小提琴的位置咱们是知道的,分别在0,1,2,3,5,6,8,9位置上,所以咱们先绘制八个癌与癌旁的小label

  1. x=c(1:5)

  2. y=c(1:5)

  3. plot(x, y, xlim=c(0,10),xlab='',ylab='Gene Expression'

  4. ,main=''

  5. , ylim=c(min(sigData),max(sigData)+1),pch=21,col='white',xaxt="n")

  6. for(i in 1:4){

  7.  basal=sigData[tumor,i]

  8.  noBasal=sigData[normal,i]

  9.  vioplot(add = T,col = 'blue',basal,at=3*(i-1),lty=1)

  10.  vioplot(add = T,col = 'red',noBasal,at=3*(i-1)+1,lty=1)

  11.  p=round(wilcox.test(basal,noBasal)$p.value,3)

  12.  mx=max(c(basal,noBasal))

  13.  lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx))

  14.  text(x=3*(i-1)+0.5,y=mx+0.5,labels=ifelse(p==0,paste0("p<0.001")

  15.  ,paste0("p,p)),cex = 0.5)

  16. }

  17. axis(side=1, at = c(0,1,3,4,6,7,9,10), labels = F, tick = TRUE)

  18. text(c(0,1,3,4,6,7,9,10),rep(par("usr")[3]-0.7,8),xpd = NA,cex = 0.8,

  19. labels=rep(c('Basal','non-Basal'),4),srt = 45,adj=1)


对比代码咱们加了两行


  1. axis(side=1, at = c(0,1,3,4,6,7,9,10), labels = F, tick = TRUE)#表示坐标抽啦,side=1表示横轴,at= 表示八个位置,labels=F表示不绘制label,tick=T表示显示那个小竖杆

  2. text(c(0,1,3,4,6,7,9,10),rep(par("usr")[3]-0.7,8),xpd = NA,cex = 0.8,labels=rep(c('Basal','non-Basal'),4),srt = 45,adj=1)#绘制那些文字啦,srt表示倾斜角度#注意:rep(par("usr")[3]-0.7,8)很重要,他直接表示了垂直方向的位置,尤其是那个0.7,可以试着改大改小,好好体会一下

最好要把基因给标上去嘛,代码如下:

  1. x=c(1:5)

  2. y=c(1:5)

  3. plot(x, y, xlim=c(0,10),xlab='',ylab='Gene Expression'

  4. ,main=''

  5. , ylim=c(min(sigData),max(sigData)+1),pch=21,col='white'

  6. ,xaxt="n")

  7. for(i in 1:4){

  8.  basal=sigData[tumor,i]

  9.  noBasal=sigData[normal,i]

  10.  vioplot(add = T,col = 'blue',basal,at=3*(i-1),lty=1)

  11.  vioplot(add = T,col = 'red',noBasal,at=3*(i-1)+1,lty=1)

  12.  p=round(wilcox.test(basal,noBasal)$p.value,3)

  13.  mx=max(c(basal,noBasal))

  14.  lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx))

  15.  text(x=3*(i-1)+0.5,y=mx+0.5,labels=ifelse(p==0,paste0("p<0.001")

  16.  ,paste0("p,p)),cex = 0.5)

  17. }

  18. axis(side=1, at = c(0,1,3,4,6,7,9,10), labels = F, tick = TRUE)

  19. text(c(0,1,3,4,6,7,9,10),rep(par("usr")[3]-0.7,8),xpd = NA,cex = 0.8

  20. ,labels=rep(c('Tumor','Normal'),4),srt = 45,adj=1)

  21. text(c(0.5,3.5,6.5,9.5),rep(par("usr")[3]-3,8),xpd = NA

  22. ,labels=colnames(sigData),cex = 0.8)