師弟：師兄，火山圖怎么得到的？我都快愁死了......
師兄：你不是學(xué)過R語言嗎？
師弟：道理倒是都懂，看書歸看書，但是實(shí)際應(yīng)用起來大腦就一片空白，完全沒有構(gòu)思......
師兄：不怕，我這個(gè)有個(gè)完整的框架，你來學(xué)習(xí)一下？

首先，什么是火山圖

火山圖是轉(zhuǎn)錄組分析中經(jīng)常能夠見到的一類統(tǒng)計(jì)圖，常用于反映差異表達(dá)基因概況。例如腫瘤組織與正常組織相比，那些基因存在過表達(dá)或被抑制，可以是編碼基因、蛋白表達(dá)水平，也可以是其它lncRNA、miRNA、circRNA等非編碼RNA分子。

火山圖本質(zhì)上是散點(diǎn)圖的一種，因其外觀形似火山（噴發(fā)），所以稱為火山圖。兩個(gè)坐標(biāo)軸中，一個(gè)代表差異倍數(shù)變化（通常為log2轉(zhuǎn)化后的fold change值），一個(gè)代表顯著性（通常為-log10轉(zhuǎn)化后的p值或p調(diào)整值），將各基因映射到圖中，并根據(jù)預(yù)先定義的閾值在圖中用不同顏色標(biāo)識(shí)顯著上、下調(diào)的基因。

文獻(xiàn)中常見火山圖舉例

差異表達(dá)基因分析

可知火山圖主要反映的是基因整體差異表達(dá)狀態(tài)。因此為了獲得火山圖，首先需要根據(jù)基因表達(dá)定量數(shù)據(jù)計(jì)算各基因在比較兩組中的差異狀況，獲得表達(dá)值的倍數(shù)變化以及顯著性信息，并定義閾值確定一些高度顯著的差異基因。

常見方法如edgeR、DESeq2等，這些都是目前識(shí)別差異表達(dá)分子的主流工具。在前面的文章中也介紹過這些工具的使用，詳情可點(diǎn)擊查看：edgeR、DESeq2。

火山圖繪制

具體的差異表達(dá)基因計(jì)算過程不是本篇的內(nèi)容，就不再細(xì)提。

總之最后獲得了類似這樣的一張統(tǒng)計(jì)表，見示例數(shù)據(jù)“control_treat.count.txt”，記錄了各基因的名稱、log2轉(zhuǎn)化后的fold change值、顯著性p值、p調(diào)整值以及標(biāo)識(shí)的顯著上、下調(diào)基因等，接下來通過該表繪制火山圖就可以了。

本文的所有測試數(shù)據(jù)和R代碼，可在文末獲取。

DESeq2的示例輸出

就R的作圖包而言，ggplot2為此提供了優(yōu)秀的作圖方案，參考以下示例。

初步繪制輪廓

首先，如上所述，火山圖就是一種特殊的散點(diǎn)圖樣式。因此可以先不要管太多的細(xì)節(jié)，先大致繪制一個(gè)輪廓來觀察下，數(shù)據(jù)是否合理，差異基因數(shù)量大致有多少等。

選擇表中的log2轉(zhuǎn)化后的fold change值（log2FoldChange），以及p調(diào)整值（padj）并進(jìn)行-log10轉(zhuǎn)化后，分別作為兩個(gè)坐標(biāo)軸，繪制散點(diǎn)圖就可以了。

對于上、下調(diào)的基因，使用不同顏色表示。

#如未安裝R包需要首先安裝
#install.packages('ggplot2')

#加載R包
library(ggplot2)

#讀取數(shù)據(jù)
dat <- read.table('control_treat.DESeq2.txt', header = TRUE, sep = '\t')

#ggplot2 作圖，橫坐標(biāo) log2FoldChange，縱坐標(biāo) -log10 轉(zhuǎn)化后的 padj
p <- ggplot(data = dat, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(aes(color = sig), size = 1) +  #繪制散點(diǎn)圖，點(diǎn)的顏色表示上下調(diào)狀態(tài)
scale_color_manual(values = c('red', 'gray', 'green4'), limits = c('up', 'none', 'down')) + #指定節(jié)點(diǎn)顏色
labs(x = 'log2 Fold Change', y = '-log10 adjust p-value', 
    title = 'control vs treat', color = '')  #x、y 軸標(biāo)題設(shè)置

p

初始樣式，先繪制一個(gè)大致輪廓，用于觀測整體數(shù)據(jù)

有關(guān)主題細(xì)節(jié)的調(diào)整

大致的輪廓出來后，如果覺得數(shù)據(jù)可觀，適合擺放到文章中，就可以進(jìn)一步的調(diào)整火山圖的細(xì)節(jié)了。

接下來就需要設(shè)法將這個(gè)火山圖變得更美觀一些，涉及了ggplot2的主題調(diào)整。ggplot2的主題樣式非常多，這里以一些簡單示例來大致展示下。

例如，標(biāo)題居中、去除網(wǎng)格線和背景、去除圖例背景、字體大小修改以及在圖中標(biāo)出差異基因篩選的閾值線等。

#theme()用于修改主題，包括修改背景色、網(wǎng)格線、字體大小等
p <- p +
theme(plot.title = element_text(hjust = 0.5, size = 14), 
    panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), 
    legend.key = element_rect(fill = 'transparent'), legend.text = element_text(size = 10), 
    axis.title = element_text(size = 12), axis.text = element_text(size = 10))

p

#xlim()和ylim()可設(shè)置坐標(biāo)展示范圍，使圖左右對稱
p <- p + xlim(-12, 12) + ylim(0, 35)

p

調(diào)整主題后的示例，比剛才美觀許多

#例如先前根據(jù) fold change ≥ 2 以及 padj < 0.01 篩選的顯著差異基因
#也就對應(yīng)了 |log2FoldChange| ≥ 1 以及 -log10(padj) > 2 的閾值
#geom_vline()和geom_hline()，在圖中繪制垂線指示閾值線
p <- p +
geom_vline(xintercept = c(-1, 1), lty = 3, color = 'black') +
geom_hline(yintercept = 2, lty = 3, color = 'black')

p

添加差異基因篩選時(shí)的閾值線

特殊基因的標(biāo)識(shí)

那么，如果有一些特別重要的基因，想在圖中標(biāo)識(shí)出它們的名稱，應(yīng)該怎樣做呢？

繼續(xù)來看這個(gè)示例，示例數(shù)據(jù)“top5_gene.txt”是選擇的5個(gè)重要基因的名稱，將根據(jù)這些名稱在圖中定位基因并將它們的名稱顯示出來。

本文的所有測試數(shù)據(jù)和R代碼，可在文末獲取。

#標(biāo)識(shí)特定的基因，使用 ggplot2 的拓展包 ggrepel 來完成
library(ggrepel)

lab <- read.table('top5_gene.txt', header = TRUE, sep = '\t')
dat_select <- subset(dat, gene_id %in% lab$gene_id)  #根據(jù)重要的基因名稱在原數(shù)據(jù)中提取它們的 log2FC 和 padj 值

p +  #并將這幾個(gè)基因名稱，根據(jù)坐標(biāo)位置添加在原圖中
geom_text_repel(data = dat_select, aes(x = log2FoldChange, y = -log10(padj), label = gene_id),
    size = 3, show.legend = FALSE)

在圖中標(biāo)識(shí)重要的基因名稱，ggrepel包的好處是可以防止標(biāo)簽重疊

師兄：如何，get到了嗎？
師弟：厲害厲害，從概念到分析再到作圖，都一并包含了，點(diǎn)贊！