什么是火山圖？
火山圖其實是一種很形象的叫法，它可以通過關注對象的落點從而直觀地展示該對象的所屬區(qū)域。其通常用于展示差異結果，比如RNA-seq差異基因展示。讀懂了“火山”火星噴射的落點橫縱坐標意義，就讀懂了火山圖：
X軸：Log2 Fold Change, Fold Change指樣本間表達量的差異倍數； 取Log2是為了讓差異較大的和差異比較小的數值在視覺上縮小距離 (e.g. 原來2倍的差異等同于1Log2FC)。一般默認取log2FC絕對值大于1為差異基因的篩選標準。
Y軸：-Log(adjust P-value), 對矯正后的P值取負對數(-log)；矯正P值為多重假設檢驗矯正過的差異顯著性P值。由于轉錄組測序的差異表達分析是對大量的基因表達值進行獨立的統(tǒng)計假設檢驗會存在假陽性問題，因此在進行差異表達分析過程中，采用了公認的Benjamini-Hochberg校正方法對原有假設檢驗得到的顯著性p值（p-value）進行校正，剔除假陽性。
紅點: p<0.05且log2FC>1的基因; 藍點：p<0.05且log2FC< -1的基因; 灰色點：FC<2,即|log2FC|<1。
R繪制火山圖
下面就進入主題，用R繪制火山圖，我們的教程小白也能看懂哦??！
R包：Enhanced Volcano Plot的實例操作解說 EnhancedVolcano by Kevin Blighe
1. 安裝
if (!requireNamespace('BiocManager', quietly = TRUE))
install.packages('BiocManager')
BiocManager::install('EnhancedVolcano')
BiocManager::install('DESeq2')
BiocManager::install("airway")
BiocManager::install("magrittr")
2. 加載數據集，先安裝R包airway，magrittr。
為了快速入門及有效展示，我們參考EnhancedVolcano的官方解說，之后會對部分參數解說，調整參數畫出長在自己審美點的火山圖。（因為R包EnhancedVolcano是更新的，如果參數有變，請參照報錯提示修改代碼）
library(EnhancedVolcano)
library(airway)
library(magrittr)
data('airway')
airway$dex %<>% relevel('untrt')
ens <- rownames(airway)
library(org.Hs.eg.db)
symbols <- mapIds(org.Hs.eg.db, keys = ens, column = c('SYMBOL'),
keytype = 'ENSEMBL')
symbols <- symbols[!is.na(symbols)]
symbols <- symbols[match(rownames(airway), names(symbols))]
rownames(airway) <- symbols
keep <- !is.na(rownames(airway))
airway <- airway[keep,]
3. DESeq2篩選差異基因
library('DESeq2')
dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds <- DESeq(dds, betaPrior=FALSE)
res <- results(dds, contrast = c('dex','trt','untrt'))
4. EnhancedVolcano繪制
volcanoplot<-EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'pvalue')
ggsave("volcanoplot.png",volcanolot)
5. 調整繪圖參數
指南中的火山圖過于凌亂，調整參數使火山圖更加簡潔
##收斂結果，會改變log2FC，不改變P指。目的是為了使火山圖
res <- lfcShrink(dds, contrast = c('dex','trt','untrt'), res=res, type = 'normal')
#去除P值為NA的基因
res<-na.omit(res)
#篩選差異顯著基因，由于該實驗中顯著且|log2FC|>1的點太多，我們制定更嚴格的篩選標準，令|log2FC|>2
celltype<-rownames(res[res$padj<0.05&abs(res$log2FoldChange)>2,])
head(celltype)
p1<-EnhancedVolcano(res,
lab = rownames(res),
selectLab=celltype,
axisLabSize=13,
title="Treat vs Untreat",
subtitle=NULL,
titleLabSize=15,
#caption = paste0('Total = ', nrow(toptable), ' variables'),
captionLabSize=11,
#drawConnectors = T,
#widthConnectors = 0.2,
#colConnectors = 'grey50',
#boxedLabels=T,
labSize=2.5,
x = 'log2FoldChange',
y = 'padj',
xlim = c(-5,5),
pointSize=2,
#shape = c(6, 6, 19, 16),
colAlpha=0.7,
gridlines.major=F,
gridlines.minor=F,
border="full",
borderWidth=1,
#boxedLabels=T,
cutoffLineType="longdash",
cutoffLineCol="red",
legendVisible=F,
pCutoff=0.05,
FCcutoff=2,
#vline=c(-2,2),
xlab = bquote(~log[2]~ 'fold change'),
ylab=bquote(~-log[10]~'adjusted p-value')
)
p2<-p1+theme(axis.text.x = element_text(color="black", size=12),
axis.text.y = element_text(color="black", size=12),
plot.title = element_text(hjust = 0.5))
ggsave(p2,file="volcanoplot2.png")
這樣的火山圖是不是錯落有致，簡潔明了，是你想象中的模樣??。
6. ggplot2繪制火山圖（使上下調顯著差異基因顯示不同的顏色）
res$threshold<-as.factor(ifelse(res$log2FoldChange >= 2,'Up',ifelse(res$padj<0.05 & res$log2FoldChange <= -2,'Down','Not')))
res<-data.frame(res)
p3<-ggplot(data=res, aes(x=log2FoldChange, y=-log10(padj), colour=threshold, fill=threshold)) +
scale_color_manual(values=c("blue", "grey","red"))+
geom_point(alpha=0.6,size=2) +
#xlim(c(-6, 6)) +
#ylim(c(0, 300)) +
theme_bw(base_size = 12, base_family = "Times") +
geom_vline(xintercept=c(-2,2),lty=1,col=c("blue","red"),lwd=0.6)+
geom_hline(yintercept = -log10(0.05),lty=2,col="black",lwd=0.6)+
theme(legend.position="right",
panel.grid=element_blank(),
legend.title = element_blank(),
legend.text= element_text(face="bold", color="black",family = "Times", size=8),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color="black", size=12),
axis.text.y = element_text(color="black", size=12),
axis.title.x = element_text(face="bold", color="black", size=12),
axis.title.y = element_text(face="bold",color="black", size=12))+
labs(x="log2 (Fold Change)",y="-log10 (p-value)",title="Treat vs Untreat")
這版的風格跟EnhancedVolcano的保持一致了，我們就是這么專一??！
如果你有什么問題，在留言區(qū)域留言哦，如果你喜歡我們的文章，請點個贊吧，還可以收藏我們的簡書賬號，從9月開始，我們會定時更新哦!
