拖延癥讓我一直疲于整理qiime2的分析過程。而且手上暫時(shí)也沒有相應(yīng)的16S數(shù)據(jù)要分析，所以打算整理整理R的畫圖腳本。而且總是百度一些畫圖細(xì)節(jié)讓我“厭倦”了，覺得自己急需整理總結(jié)。本來想從最基本的條形圖、箱型圖開始，但是今天恰好做了火山圖，那我們就一起去看“火山噴發(fā)”吧。
什么是火山圖
火山圖，是形如火山噴發(fā)的一種圖形展示方法，常被用于展示差異，比如差異基因、差異微生物等等。
上圖就是一張典型的火山圖，描述了差異基因的情況。該火山圖的y軸是-log10(Qvalue)，即qvalue（pvalue校正后的值）取-log10，因此數(shù)值越高說明qvalue越小即越顯著。橫坐標(biāo)是Log2 fold change，即對(duì)fold change取log2，所以越靠?jī)蓚?cè)的點(diǎn)（每個(gè)點(diǎn)代表一個(gè)基因），其基因表達(dá)量上調(diào)或者下調(diào)幅度越大。
一般來說在差異基因分析過程中，我們通常認(rèn)為qvalue小于0.05且foldchange的絕對(duì)值大于2為差異基因。當(dāng)然在差異基因數(shù)量過多的時(shí)候，我們可以調(diào)整篩選標(biāo)準(zhǔn)以獲得相對(duì)適合數(shù)量的差異基因再進(jìn)行下游的富集分析，比如要求qvalue小于0.001，或者要求foldchange的絕對(duì)值大于5等等。當(dāng)差異基因數(shù)量過少的時(shí)候，我們可以考慮將foldchange的絕對(duì)值變?yōu)?.5，或者考慮選擇pvalue小于0.05。
而圖中的虛線就是根據(jù)自己的篩選標(biāo)準(zhǔn)確定添加。其中兩條豎線（x=-2和x=2）說明該篩選標(biāo)準(zhǔn)是要求foldchange的絕對(duì)值大于4。橫線（大膽猜測(cè)是在y=2處），說明要求qvalue小于0.01。
當(dāng)Foldchange沒有那么明顯的時(shí)候，我們的橫軸也可以選擇展示Foldchange，不取log2。
火山圖要怎么畫
(1) 需要什么格式的數(shù)據(jù)
從網(wǎng)上找了一個(gè)別人的RNA-seq數(shù)據(jù) ，讓我們一起來瞧瞧。
數(shù)據(jù)部分
第一列是基因名，第二列是計(jì)算好的log2FoldChange（通常我們用DESeq等軟件分析得到的結(jié)果里包含這一個(gè)值），第三列是pvalue，第四列是校正后的p即padj。
如果你已經(jīng)確定要畫pvalue或者padj那只要三列即可。
（2）如何使用ggpot2做火山圖
能夠做火山圖的方法有很多，有一些RNA-seq分析的包中自帶了畫火山圖的函數(shù)。利用R自帶的基礎(chǔ)畫圖函數(shù)也可以畫，但是鑒于之后我們都幾乎都選擇ggplot2包進(jìn)行作圖，所以只展示如何用ggplot2包畫圖。
#加載包
library(ggplot2)
#讀取數(shù)據(jù)
Dat<-read.table('./results.txt',header = T,stringsAsFactors = F,check.names = F,sep=' ')
#作圖
ggplot(Dat,aes(x=log2FoldChange,y=-log10(padj)))+
geom_point()
通過簡(jiǎn)單的幾行代碼我們就能夠得到一個(gè)最基本（簡(jiǎn)陋）的火山圖。
然后我們來美化這張火山圖。
#加載包
library(ggplot2)
library(ggrepel)
#讀取數(shù)據(jù)
Dat<-read.table('./results.txt',header = T,stringsAsFactors = F,check.names = F,sep=' ')
#確定是上調(diào)還是下調(diào)，用于給圖中點(diǎn)上色）
Dat$threshold = factor(ifelse(Dat$padj < 0.05 & abs(Dat$log2FoldChange) >= 1, ifelse(Dat$log2FoldChange>= 1 ,'Up','Down'),'NoSignifi'),levels=c('Up','Down','NoSignifi'))
ggplot(Dat,aes(x=log2FoldChange,y=-log10(padj),color=threshold))+
geom_point()+
scale_color_manual(values=c("#DC143C","#00008B","#808080"))+#確定點(diǎn)的顏色
geom_text_repel(
data = Dat[Dat$padj<0.05&abs(Dat$log2FoldChange)>1,],
aes(label = Gene),
size = 3,
segment.color = "black", show.legend = FALSE )+#添加關(guān)注的點(diǎn)的基因名
theme_bw()+#修改圖片背景
theme(
legend.title = element_blank()#不顯示圖例標(biāo)題
)+
ylab('-log10 (p-adj)')+#修改y軸名稱
xlab('log2 (FoldChange)')+#修改x軸名稱
geom_vline(xintercept=c(-1,1),lty=3,col="black",lwd=0.5) +#添加橫線|FoldChange|>2
geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=0.5)#添加豎線padj<0.05
Volcano
一張精致的火山圖就做好啦。
其中說明一下lty是用于描述添加線的類型，數(shù)字對(duì)應(yīng)不同的類型，具體的可以參考下圖：
