什么是火山圖 

  火山（Volcano Plot）圖在一張圖中顯示了兩個(gè)重要的指標(biāo)（Foldchange/pvalue），可以非常直觀且合理地篩選出在兩樣本間發(fā)生差異表達(dá)的基因。檢驗(yàn)分析出兩樣本間顯著差異表達(dá)的基因后，以log2（fold change）為橫坐標(biāo)，以T檢驗(yàn)顯著性檢驗(yàn)P值的負(fù)對(duì)數(shù)-log10（pvalue）為縱坐標(biāo)，即可得火山圖（Volcano Plot）。

火山圖怎么制作呢  

火山圖（Volcano Plot）可以直觀地展現(xiàn)所有基因的FDR與Fold Change之間的關(guān)系，以便快速查看基因在兩組樣品間的表達(dá)水平差異程度及其統(tǒng)計(jì)學(xué)顯著性。

兩個(gè)重要概念：FDR和FC

FDR：在差異表達(dá)分析過程中，采用Benjamini-Hochberg方法對(duì)原有假設(shè)檢驗(yàn)得到的顯著性p值（p-value）進(jìn)行校正，并最終采用校正后的p值，即FDR（False Discovery Rate）。

FC（Fold Change）：兩樣品（組）間表達(dá)量的比值。

我們簡(jiǎn)單演示一下

#載入包

library(ggplot2)

#導(dǎo)入數(shù)據(jù)

df<-read.delim(file = 'test.xls',header = T,sep = '\t',check.names = F,row.names = 1,stringsAsFactors = F)

#簡(jiǎn)單查看

head(df,n=10)

針對(duì)數(shù)據(jù)進(jìn)行處理，在這里根據(jù)兩（組）樣品之間表達(dá)水平的相對(duì)高低，差異表達(dá)基因可分為上調(diào)基因（upregulated Gene）和下調(diào)基因（down）。上調(diào)和下調(diào)是相對(duì)的，由所給A和B的順序決定，若更換A和B的順序會(huì)上調(diào)基因和下調(diào)基因反過來，但對(duì)結(jié)果沒有影響。

我們根據(jù)FDR<=0.05和|FC|>=2作為篩選標(biāo)準(zhǔn)！

for (i in 1:nrow(df) ) {
  if ( as.double(df[i,ncol(df)-1]) >= 2 && as.double(df[i,ncol(df)-2]) <= 0.05 ) {
    df[i,ncol(df)] <- 'up'
  } else if ( as.double(df[i,ncol(df)-1]) <= -2 && as.double(df[i,ncol(df)-2]) <= 0.05 ) {
    df[i,ncol(df)] <- 'down'
  } else {
    df[i,ncol(df)] <- 'normal'
  }
}

然后呢，就是利用ggplot2進(jìn)行繪制了……

p <- ggplot(df, aes(x=log2FC, y=-log10(FDR), colour=regulated) ) geom_point(size=1)
p <- p scale_color_manual(values = c('red','green','blue'))
p <- p geom_vline(xintercept=c(-2,2), linetype='longdash', size=0.2)
p <- p geom_hline(yintercept=c(-log10(0.05)), linetype='longdash', size=0.2)
p

差異表達(dá)火山圖中的每一個(gè)點(diǎn)表示一個(gè)基因，橫坐標(biāo)表示某一個(gè)基因在兩樣品中表達(dá)量差異倍數(shù)的對(duì)數(shù)值，其絕對(duì)值越大，說明表達(dá)量在兩樣品間的表達(dá)量倍數(shù)差異越大；縱坐標(biāo)表示FDR的負(fù)對(duì)數(shù)值，其值越大，表明差異表達(dá)越顯著，篩選得到的差異表達(dá)基因越可靠。圖中紅色和藍(lán)色的點(diǎn)代表有顯著性表達(dá)差異的基因，紅色代表基因表達(dá)量下調(diào)，藍(lán)色代表基因表達(dá)量上調(diào)，綠色的點(diǎn)代表無顯著性表達(dá)差異的基因。

拓展

MA圖可以直觀地查看兩組樣品中基因的表達(dá)豐度和差異倍數(shù)的整體分布！

#統(tǒng)計(jì)上調(diào)基因和下調(diào)基因的數(shù)目

normal<-paste0('normal:',sum(df[,'regulated'] == 'normal'))
up<-paste0('up:',sum(df[,'regulated'] == 'up'))
down<-paste0('down:',sum(df[,'regulated'] == 'down'))
sig <- c(normal,up,down)

#繪制MA圖

v<-c('sample1','sample2')
data=data.frame(x=log2(rowMeans(df[,v])), y=df[,c('log2FC')],lab=factor(df[,'regulated'],levels=unique(df[,'regulated'])))
p=ggplot(data,mapping=aes(x=log2(rowMeans(df[,v])), y=df[,c('log2FC')],color=lab))
p=p geom_point(size=1)
p=p xlab('log2(FPKM)') ylab('log2(FC)') ggtitle('MA plot') theme(plot.title=element_text(face='bold',size=14))
p=p theme_classic()
p=p scale_color_manual(name='regulated',values=c('red','green','blue'),labels=sig)
p

這里我們提供另外一種繪制火山圖的方法：

library(ggplot2)

data=read.table(file='deg-result.txt',header=T,row.names=1)

threshold<-as.factor((datalogFC>1.5|data

logFC< -1.5)& data$P.Value<0.05)

ggplot(data,aes(x=logFC,y=-log10(P.Value),colour=threshold)) xlab('log2fold change') ylab('-log10p-value') geom_point()

關(guān)于p值