【科研貓·繪圖】朋友圈最火熱的「火山圖」（帶R程式碼分享）

2019 年 10 月 6 日
筆記

暑假到了，在這個火熱的夏日，還有比一張美麗的文章主圖更能讓我們充滿激情和力量嗎？在這組學撐起生物資訊學半邊天的時代，火山圖，在轉錄組、蛋白組、代謝組等多種組學數據中，都佔有舉足輕重的地位。

那麼，什麼是火山圖（Volcano Plot）。正如名稱所描述，火山圖的圖形非常像火山噴發的形狀。火山圖通常用來展示差異表達的基因，常常出現在晶片、轉錄組、蛋白組、代謝組等組學檢測技術的結果中，並且通常伴隨熱圖一起出現。

火山圖是散點圖的一種展現形式。以實驗組和對照組轉錄組數據為例，一張完美的火山圖，通常由幾個部分組成，顯著上調差異表達基因，顯著下調差異表達基因。一般來說，x軸為實驗組基因表達量比上對照組基因表達量的倍數差異，而y軸則為實驗組比對照組之後的p值或者校正後的p值。火山圖上，一個點代表一個基因，而顏色則代表他們是顯著上調還是顯著下調。

好了，基礎知識就介紹到這裡，下面讓我們來看如何繪製一張完美的火山圖。

首先，讓我們先看一下數據。我們的數據使用的是RNAseq表達譜計算出來的差異基因，差異基因的計算方法為limma。在RNAseq分析中，使用較多的計算差異基因的軟體為DESeq2和limma。差異分析結果文件我們存儲在DEGdata.txt文件中，使用excel打開改文件，一共有7列數據：第一列Symbol為基因名；第二列logFC為實驗組比上對照組的倍數差異，其中大於0為實驗組上調基因，小於0為實驗組下調基因；第三列AveExpr為該基因在所有樣本中表達量均值；第四列t為t統計量；第五列P.value為實驗組比對照組的P值；第六列adj.p.value為校正後的P值，因為基因和基因並不是相互獨立的，所以我們需要對P值進行校正來降低結果的假陽性，常用的校正方法為FDR校正；第七列B為B統計量。

然後，讓我們打開Rstudio，開始安裝一些必須的包，這裡推薦一款ggplot2一鍵繪圖神器，ggpubr (http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/)。這個包是在ggplot2的基礎上，整合了很多繪圖函數，對用戶十分友好。另外需要用到的一個包圍ggthemes，這個包整合了很多ggplot2的繪圖主題，對科研人員十分友好。

在包安裝完成之後，我們將它們載入到R環境中。接下來就是讀取我們差異基因的文件了。差異基因文件我們存儲在DEGdata.txt文件中。

在R中查看文件前6行。

這7列文件之前已經介紹過，我們畫火山圖，只需要其中的logFC和adj.P.val就可以了。在繪圖之前，我們需要對adj.P.val進行轉換，將它的值變成-1 * log10，這樣的話可以拉開差異表達基因之間的間距。

接下來就是繪製火山圖了~ 程式碼如下：

畫出來的火山圖是這樣色兒的：

當然，直接這麼畫火山圖，結果是非常丑的。編輯們看到了，估計會直接拒稿吧(╥╯^╰╥)，所以我們需要對它進行美化~

首先需要區分其中哪些是顯著差異表達基因。因此，我們需要對logFC和adj.P.val兩列進行過濾。在本文中，我們設置的過濾的條件為，adj.P.val小於0.05並且logFC大於2（4倍差異）為顯著上調差異表達基因，adj.P.val小於0.05並且logFC小於-2（4倍差異）為顯著下調差異表達基因。大家可以根據實驗結果，適當調整logFC。