为了表示基因在样本中的差异,对于许多个基因的话用火山图感觉挺高大上的,那么对于少数几个基因的话一般用什么来展示其表达差异呢,有很多种办法,比如箱线图,柱状图等等,今天主要用小提琴图来可视化几个基因的差异表达。
首先安装R包
source("https://bioconductor.org/biocLite.R")
biocLite('vioplot')
然后我的数据是一个矩阵,行为样本列为基因如:
已经知道癌症样本和正常样本的索引
那么绘制小提琴图的简单办法如下:
vioplot(sigData[tumor,1]
,sigData[normal,1]
,sigData[tumor,2]
,sigData[normal,2]
,sigData[tumor,3]
,sigData[normal,3]
,sigData[tumor,4]
,sigData[normal,4]
,names=rep(c("Tumor","Normal"),4),col="red")
从图中可以看出这四个基因其实在癌与癌旁是有差异的,但是好丑,很显然横轴要标上基因,每一个基因的癌与癌旁近一点,同时标记上p值最好了
那么怎么调整呢
首先需要调整的小提琴图的间距,让其看起来像四组,而不是八组,我们机智的设计了x轴,如代码所示:
x=c(1:5)
y=c(1:5)
plot(x, y, xlim=c(0,10),xlab='',ylab='Gene Expression'
,main=''
, ylim=c(min(sigData),max(sigData)+1),pch=21,col='white',xaxt="n")
for(i in 1:4){
basal=sigData[tumor,i]
noBasal=sigData[normal,i]
vioplot(add = T,basal,at=3*(i-1),lty=1)
vioplot(add = T,noBasal,at=3*(i-1)+1,lty=1)
}
从代码中可以看出我们将x轴设计成了0-10,然后再0,1,2,3,5,6,8,9处画上小提琴,这样就可以看出四组数据了,顺便将y轴的label加上
但是这个还是很丑,我们想将tumor和normal用颜色区分一下下,怎么办,如代码:
x=c(1:5)
y=c(1:5)
plot(x, y, xlim=c(0,10),xlab='',ylab='Gene Expression'
,main=''
, ylim=c(min(sigData),max(sigData)+1),pch=21,col='white',xaxt="n")
for(i in 1:4){
basal=sigData[tumor,i]
noBasal=sigData[normal,i]
vioplot(add = T,col = 'blue',basal,at=3*(i-1),lty=1)
vioplot(add = T,col = 'red',noBasal,at=3*(i-1)+1,lty=1)
}
如咱们所愿,颜色变了,但是p值木有加上,还是挺遗憾的,咱们使用‘Mann-Whitney’ test来检验差异显著性,然后把p值在图中,做法如下:
x=c(1:5)
y=c(1:5)
plot(x, y, xlim=c(0,10),xlab='',ylab='Gene Expression'
,main=''
, ylim=c(min(sigData),max(sigData)+1),pch=21,col='white',xaxt="n")
for(i in 1:4){
basal=sigData[tumor,i]
noBasal=sigData[normal,i]
vioplot(add = T,col = 'blue',basal,at=3*(i-1),lty=1)
vioplot(add = T,col = 'red',noBasal,at=3*(i-1)+1,lty=1)
p=round(wilcox.test(basal,noBasal)$p.value,3)
mx=max(c(basal,noBasal))
lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx))
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值已经加上去了,对比之前的代码主要多了四行代码
p=round(wilcox.test(basal,noBasal)$p.value,3)#用来检验基因的差异显著性p值,并去3位有效数字
mx=max(c(basal,noBasal))#找到新画上去的两个小提琴的最高的位置,为画上p值的位置坐准备
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
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值的文本啦,x,y都是图中位置,仔细体会,label就是文本,cex设成0.5表示原来的50%
p值也画上去了,但是x轴一直没有,最后咱们来绘制x轴,首先分析一下八个小提琴的位置咱们是知道的,分别在0,1,2,3,5,6,8,9位置上,所以咱们先绘制八个癌与癌旁的小label
x=c(1:5)
y=c(1:5)
plot(x, y, xlim=c(0,10),xlab='',ylab='Gene Expression'
,main=''
, ylim=c(min(sigData),max(sigData)+1),pch=21,col='white',xaxt="n")
for(i in 1:4){
basal=sigData[tumor,i]
noBasal=sigData[normal,i]
vioplot(add = T,col = 'blue',basal,at=3*(i-1),lty=1)
vioplot(add = T,col = 'red',noBasal,at=3*(i-1)+1,lty=1)
p=round(wilcox.test(basal,noBasal)$p.value,3)
mx=max(c(basal,noBasal))
lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx))
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)
}
axis(side=1, at = c(0,1,3,4,6,7,9,10), labels = F, tick = TRUE)
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)
对比代码咱们加了两行
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表示显示那个小竖杆
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,可以试着改大改小,好好体会一下
最好要把基因给标上去嘛,代码如下:
x=c(1:5)
y=c(1:5)
plot(x, y, xlim=c(0,10),xlab='',ylab='Gene Expression'
,main=''
, ylim=c(min(sigData),max(sigData)+1),pch=21,col='white'
,xaxt="n")
for(i in 1:4){
basal=sigData[tumor,i]
noBasal=sigData[normal,i]
vioplot(add = T,col = 'blue',basal,at=3*(i-1),lty=1)
vioplot(add = T,col = 'red',noBasal,at=3*(i-1)+1,lty=1)
p=round(wilcox.test(basal,noBasal)$p.value,3)
mx=max(c(basal,noBasal))
lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx))
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)
}
axis(side=1, at = c(0,1,3,4,6,7,9,10), labels = F, tick = TRUE)
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('Tumor','Normal'),4),srt = 45,adj=1)
text(c(0.5,3.5,6.5,9.5),rep(par("usr")[3]-3,8),xpd = NA
,labels=colnames(sigData),cex = 0.8)