除了火山圖,差異表達基因還可以這樣展示

R語言繪製基因表達基因的「對稱散點圖」

 

轉錄組分析中，計算了兩組間差異表達的基因後，通常怎樣表示？您可能第一時間想到可以使用火山圖。的確，火山圖是使用頻率最多的，在火山圖中可以很輕鬆地根據基因在兩組間的Fold Change值以及顯著性p值，識別和判斷差異表達基因概況。火山圖實質上就是一種散點圖，通常橫縱坐標分別代表了log2轉化後的Fold Change以及-log10轉化後的p值或p調整值信息（下圖左）。提到散點圖，常見的還有另一種展示差異表達基因的樣式：橫縱坐標軸可分別代表兩組基因表達均值，這種風格可以更方便直觀對比基因在兩組中的差異狀態。

本篇教程就讓我們來學習如何繪製右圖這種「對稱散點圖」，展示組間差異基因表達格局。示例數據和R代碼等，可點擊下方「閱讀原文」獲取。（備註：若無法成功打開連結，請切換到電腦端網頁點擊下載）

示例文件「gene_diff.txt」是一組基因差異表達分析結果，記錄了處理組（treat）和對照組（control）間表達顯著不一致的基因，鑑定標準為p<0.01以及|log2 Fold Change|≥1。

其中，gene_id為基因名稱；control和treat代表了兩組中基因的平均表達值；log2FoldChange即log2轉化後的基因表達差異倍數；pvalue是差異基因顯著性p值；diff為根據p<0.01以及|log2 Fold Change|≥1篩選的差異基因，該列中「up」為上調，「down」為下調，「none」為非差異基因。

接下來通過該示例文件，展示使用R語言繪製差異基因表達「對稱散點圖」過程。

首先對數據做一些預處理。

例如，基因表達值數量級相差過大，取個對數轉換；基因名稱按是否為差異基因作個排序，避免後續作圖時被不顯著的基因點遮蓋，即排序的目的是讓這些顯著基因的點都位於圖的上方。

express <- read.delim('gene_diff.txt', sep = '\t')
express$control <- log(express$control+1)express$treat <- log(express$treat+1)
express$diff <- factor(express$diff, levels = c('up', 'down', 'none'))express <- express[order(express$diff, decreasing = TRUE), ]
head(express)  

下來就可以使用預處理後的數據作圖了。
第一種類型是將基因按上調、下調或不顯著類型著色，便於從圖中辨認差異基因。我們使用ggplot2的方法繪製差異基因散點圖。
library(ggplot2)
ggplot(express, aes(x = control, y = treat)) +geom_point(aes(color = diff), size = 1) +  scale_color_manual(values = c('red', 'gray', 'green4'), limit = c('up', 'none', 'down')) +  theme_bw() +  labs(x = 'control group', y = 'treat group', color = '') +  geom_abline(intercept = 1, slope = 1, col = 'black', linetype = 'dashed', size = 0.5) +  geom_abline(intercept = -1, slope = 1, col = 'black', linetype = 'dashed', size = 0.5) +geom_abline(intercept = 0, slope = 1, col = 'black', linetype = 'dashed', size = 0.5)
兩個坐標軸分別代表了處理組（treat）和對照組（control），圖中的點代表各基因在兩組中的平均表達值（已經作了log轉換）。treat組和control組相比，上調基因以紅色表示，下調基因以綠色表示。圖中的虛線代表了|log2FC|=1時的閾值線。
在該圖中，我們可以很輕鬆地觀察差異基因整體分布狀態和數量比較的信息。

上圖中沒有將p值信息展示出。因此另一種思路是，顏色代表p值，這樣就可以在圖中獲得一個漸變梯度。同樣使用ggplot2的方法繪製，和上述過程相比僅在顏色指定上存在區別。
ggplot(express, aes(x = control, y = treat)) +geom_point(aes(color = pvalue), size = 0.8) +  scale_color_gradient2(low = 'red', mid = 'darkgoldenrod2', high = 'royalblue2', midpoint = 0.5) +  theme_bw() +  labs(x = 'control group', y = 'treat group', color = 'p-value') +  geom_abline(intercept = 1, slope = 1, col = 'black', linetype = 'dashed', size = 0.5) +  geom_abline(intercept = -1, slope = 1, col = 'black', linetype = 'dashed', size = 0.5) +geom_abline(intercept = 0, slope = 1, col = 'black', linetype = 'dashed', size = 0.5)
類似上圖，兩個坐標軸分別代表了處理組（treat）和對照組（control），圖中的點代表各基因在兩組中的平均表達值（已經作了log轉換），圖中的虛線代表了|log2FC|=1時的閾值線。
和上圖不同點在於，此時基因按顯著性p值著色，從不顯著>顯著展示以藍色>紅色漸變，就獲得了一種梯度信息。這樣可以很方便地看出，在兩組中的表達值差異越大的基因，p值越小，二者趨勢是一致的，重在描述了差異倍數和p值的關係。
此外，若老師或同學們有RNAseq（mRNA、lncRNA、miRNA、circRNA）或蛋白質組等數據分析、繪圖等問題疑問，歡迎掃描下方二維碼回復，我們會根據大家的需求，選擇合適的問題，整理教程。
上海生因生物有著豐富的轉錄組測序、外顯子測序數據分析的經驗，同時還提供文獻或分析思路整理、GEO、TCGA公共數據挖掘、高級個性化定製分析等服務。有這方面試驗或數據分析需要的老師，可以添加技術微信聯繫我們，共同探討如何尋找基因、分子研究，如何確定分子機制。對於已經在我們公司做過測序的老師，或者打算即將在我們公司做測序的老師，可以享受免費的售後分析服務。
關注技術微信聯繫數據分析
長按識別二維碼諮詢實驗、分析
李紀偉丨寫
趙青海丨審
其他相關資料
視頻講解：文獻中常見的信號通路是如何富集出來的
視頻講解：沒想到基因GO富集分析這麼簡單
視頻講解-GEO測序數據下載軟體prefech使用視頻（windows系統）
視頻講解-veen圖在線繪製教程
R語言作圖-蜂群圖，讓組間基因表達值的比較更優雅
R語言作圖-R語言繪製「密度提琴圖」，讓提琴圖更加豐富多彩
R語言作圖-R語言繪製ceRNA網絡的衝擊圖（桑基圖）
R語言作圖-R語言繪製基因表達相關性弦圖教程
R語言作圖-R語言繪製基因表達熱圖並設置基因展示範圍
R語言作圖-R語言繪製配對箱線圖
點擊閱讀原文查看更多信息