R語言繪圖 | R繪製火山圖 EnhancedVolcano+ggplot

什麼是火山圖？

火山圖其實是一種很形象的叫法，它可以通過關注對象的落點從而直觀地展示該對象的所屬區域。其通常用於展示差異結果，比如RNA-seq差異基因展示。讀懂了「火山」火星噴射的落點橫縱坐標意義，就讀懂了火山圖：

圖片1.png

X軸：Log2 Fold Change, Fold Change指樣本間表達量的差異倍數；取Log2是為了讓差異較大的和差異比較小的數值在視覺上縮小距離  (e.g. 原來2倍的差異等同於1Log2FC)。一般默認取log2FC絕對值大於1為差異基因的篩選標準。

Y軸：-Log(adjust P-value), 對矯正後的P值取負對數(-log)；矯正P值為多重假設檢驗矯正過的差異顯著性P值。由於轉錄組測序的差異表達分析是對大量的基因表達值進行獨立的統計假設檢驗會存在假陽性問題，因此在進行差異表達分析過程中，採用了公認的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("magrittr")  BiocManager::install('EnhancedVolcano')    BiocManager::install('DESeq2')    BiocManager::install("airway")
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) <- symbolskeep <- !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的基因

#篩選差異顯著基因，由於該實驗中顯著且|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")


來源：https://www.jianshu.com/p/e2a97bdee27e
BioMan主要報導生命科學領域熱點資訊、解讀前沿進展、分享科研資料。我們組建了10餘個交流群，歡迎大家進群交流。添加公眾號博主微信：mBioMan（下方二維碼），邀你進群。溫馨提示：添加博主時，請備註一下研究方向+單位/學校！