師弟：師兄，火山圖怎么得到的格侯？我都快愁死了......
師兄：你不是學(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)識顯著上近范、下調(diào)的基因。
差異表達(dá)基因分析
可知火山圖主要反映的是基因整體差異表達(dá)狀態(tài)延蟹。因此為了獲得火山圖评矩，首先需要根據(jù)基因表達(dá)定量數(shù)據(jù)計(jì)算各基因在比較兩組中的差異狀況，獲得表達(dá)值的倍數(shù)變化以及顯著性信息阱飘，并定義閾值確定一些高度顯著的差異基因斥杜。
常見方法如edgeR、DESeq2等沥匈，這些都是目前識別差異表達(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)識的顯著上师崎、下調(diào)基因等默终，接下來通過該表繪制火山圖就可以了。
本文的所有測試數(shù)據(jù)和R代碼，可在文末獲取齐蔽。
就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
有關(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
#例如先前根據(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
特殊基因的標(biāo)識
那么盈厘，如果有一些特別重要的基因睁枕，想在圖中標(biāo)識出它們的名稱，應(yīng)該怎樣做呢沸手？
繼續(xù)來看這個(gè)示例外遇，示例數(shù)據(jù)“top5_gene.txt”是選擇的5個(gè)重要基因的名稱，將根據(jù)這些名稱在圖中定位基因并將它們的名稱顯示出來契吉。
本文的所有測試數(shù)據(jù)和R代碼跳仿，可在文末獲取。
#標(biāo)識特定的基因栅隐，使用 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)
師兄：如何租悄，get到了嗎？
師弟：厲害厲害恩袱，從概念到分析再到作圖泣棋，都一并包含了，點(diǎn)贊畔塔！
