整理思路:
先畫(huà)出中間的火山圖,再畫(huà)左右兩張放大版的圖毙玻,最后用AI將三張圖拼接到一起勺届,鋼筆工具畫(huà)出中間連接的小T形,同色的小矩形把連接處的實(shí)線(xiàn)蓋住枣购,打造出一種藕斷絲連的神秘感(偽嬉探。
1. 數(shù)據(jù)輸入和清洗
數(shù)據(jù)的輸入為 DESeq2的輸出 'res'
接下來(lái)進(jìn)行一些清洗,整理成畫(huà)圖的形式
#增加基因名一列
res<-cbind(rownames(res),res)
colnames(res)[1]<-'Symbol'
#增加不同的顏色棉圈,將P矯正<0.05且log2FC的絕對(duì)值>1作為閾值涩堤,表達(dá)上調(diào)為紅色,表達(dá)下調(diào)為藍(lán)色
a <- dim(res)[1]
res$group<-sapply(1:a,function(i){if(res$log2FoldChange[i]>1 & res$padj[i] < 0.05){"#dd6a64"}
else if(res$log2FoldChange[i]<(-1) & res$padj[i] < 0.05){"#6d95e6"}else{"grey40"}})
#增加不同點(diǎn)的尺寸分瘾,這里可以根據(jù)自己的喜好自定義
res$size<-sapply(1:a,function(i){if(abs(res$log2FoldChange[i])>1 & res$padj[i] < 0.05){"2"}
else if(abs(res$log2FoldChange[i])>1 & abs(res$log2FoldChange[i])<1.5 & res$padj[i] < 0.05){"2.5"}
else if(abs(res$log2FoldChange[i])>1.5 & abs(res$log2FoldChange[i])<2 & res$padj[i] < 0.05){"3"}
else if(abs(res$log2FoldChange[i])>2 & res$padj[i] < 0.05){"3.5"}
else{"1.5"}})
2. 畫(huà)中間的plot
fc = res$log2FoldChange
p = -log10(res$padj) #可以將P矯正換成P胎围,取決于你實(shí)驗(yàn)的設(shè)計(jì)
sizes = as.numeric(res$size)
pdf('test1.pdf')
plot(x= fc,y= p,log='y', # indicating if x or y or both coordinates should be plotted in log scale.
pch=16,
xlab=bquote(~Log[2]~'Foldchange'),
ylab=bquote(~Log[10]~'Padj'),
cex=sizes,
col=res$group,
xlim=range(fc*1.2))
至此主圖的框架已畫(huà)完,下面是添油加醋,讓他變得更好看
###
## 增加水平線(xiàn)和垂直線(xiàn)
abline(h=-log10(0.05), lty=2, lwd=1)
abline(h=2, lty=3, lwd=1)
abline(v=-1, col="blue", lty=2, lwd=1)
abline(v=1, col="red", lty=2, lwd=1)
## 增加點(diǎn)外面的黑圈
m<-which(p > -log10(0.05) & abs(fc)>1)
points(x=fc[m],y=p[m],pch=1,cex=sizes[m])
## 在主圖左右兩側(cè)畫(huà)一個(gè)黃色的長(zhǎng)方形
# 先要獲得橙色的RGB參數(shù)
rgb_orange<-col2rgb('orange')/255
# 獲得透明度為0.1時(shí)的顏色代碼
color_rect<-rgb(rgb_orange[,1][1],rgb_orange[,1][2],rgb_orange[,1][3],alpha=0.1)
# 加上長(zhǎng)方形
# rect(xleft, ybottom, xright, ytop,XXX)
rect(-3.064454,-log10(0.05),-1,12, col = color_rect)
rect(1,-log10(0.05),2.915735,12, col = color_rect)
dev.off()
注1:在加長(zhǎng)方形時(shí)白魂,四個(gè)角的坐標(biāo)不太容易獲得汽纤,可以往下到最后看我的方法
3. 畫(huà)左右兩個(gè)小圖
先將上調(diào)表達(dá)和下調(diào)表達(dá)的挑出來(lái),分別畫(huà)
res1<-res %>% dplyr::filter(log2FoldChange < (-1) & padj < 0.05)
pdf('test2.pdf',width=5,height=5)
plot(res1$log2FoldChange, -log10(res1$padj), log='y',
col=res1$group,
pch=16,
ylab=bquote(~-Log[10]~"Padj value"), xlab=bquote(~Log[2]~'Foldchange'),
cex=as.numeric(res1$size),
xlim=c(-2.5,-1),
ylim=c(1.5,6))
points(res1$log2FoldChange, -log10(res1$padj), pch=1,cex=as.numeric(res1$size))
text(res1$log2FoldChange, -log10(res1$padj), res1[,1],pos=2, #1, 2, 3 and 4, respectively indicate positions below, to the left of, above and to the right of the specified coordinates.
col='#6d95e6',cex = 1)
rect(-2.7,1,-0.7,6.4, col = color_rect)
dev.off()
有的基因名重疊了碧聪,這里可以將離得比較近的點(diǎn)挑出來(lái)冒版,換一個(gè)方向加基因名,但是我懶就不寫(xiě)了
同理逞姿,右圖也以此類(lèi)推
res2<-res %>% dplyr::filter(log2FoldChange > 1 & padj < 0.05)
pdf('red_heat.pdf',width=5,height=5)
#cols.alpha <- add.alpha(res[,7], alpha=0.6)
plot(res2$log2FoldChange, -log10(res2$padj), log='y',
col=res2$group,
pch=16, #實(shí)心圓點(diǎn)
ylab=bquote(~-Log[10]~"Padj value"), xlab=bquote(~Log[2]~'Foldchange'),
cex=as.numeric(res2$size), # 用小泡泡畫(huà)不感興趣的pathway
xlim=c(1,2.5))
points(res2$log2FoldChange, -log10(res2$padj), pch=1,cex=as.numeric(res2$size))
text(res2$log2FoldChange, -log10(res2$padj), res2[,1],pos=4, #1, 2, 3 and 4, respectively indicate positions below, to the left of, above and to the right of the specified coordinates.
col='#dd6a64',cex = 1)
gray.alpha <- add.alpha("orange", alpha=0.1)
rect(0.7,1,3,10, col = gray.alpha)
dev.off()
最后用AI將三個(gè)圖拼接辞嗡,鋼筆工具很簡(jiǎn)單就能畫(huà)出中間連接的效果,你會(huì)比我畫(huà)的更好滞造,我懶得做細(xì)活续室,草草畫(huà)了一下,還能看到毛邊谒养。
此文版權(quán)歸作者所有挺狰,未經(jīng)允許嚴(yán)禁轉(zhuǎn)載。特此警告某某某买窟,某某某丰泊,和某某某。
(算了我就是說(shuō)說(shuō)始绍,我知道你們還是會(huì)未經(jīng)允許轉(zhuǎn)載的瞳购,佛了佛了……東京今天的陽(yáng)光真好,還有什么可值得計(jì)較的呢亏推?快出去走走吧学赛!
注:rect函數(shù)中四個(gè)角的參數(shù),可用完plot函數(shù)后吞杭,用par('usr')盏浇,就會(huì)返回4個(gè)角的坐標(biāo)。但是不是很準(zhǔn)芽狗,因?yàn)閜lot的坐標(biāo)不是從整點(diǎn)開(kāi)始的绢掰,就需要配合自己的嘗試。另外童擎,在plot函數(shù)中曼月,我們添加了log='y'讓圖形變得更好看,但這樣就不會(huì)準(zhǔn)確的返回坐標(biāo)值了柔昼,所以建議重新跑一下沒(méi)有l(wèi)og='y'的plot哑芹,再par('usr')。說(shuō)的有點(diǎn)亂捕透,自己試一下就知道啦聪姿!