R語言學習 - 散點圖繪制

散點圖

散點圖在生物信息分析中是應(yīng)用比較廣的一個圖，常見的差異基因火山圖、功能富集分析泡泡圖、相關(guān)性分析散點圖、抖動圖、PCA樣品分類圖(后續(xù)推出)等。凡是想展示分布狀態(tài)的都可以用散點圖。

橫縱軸都為數(shù)字的散點圖解析

繪制散點圖的輸入一般都是規(guī)規(guī)矩矩的矩陣，可以讓不同的列分別代表X軸、Y軸、點的大小、顏色、形狀、名稱等。

輸入數(shù)據(jù)格式 (使用火山圖的輸入數(shù)據(jù)為例)

火山圖需要的數(shù)據(jù)格式如下

• id: 不是必須的，但一般的軟件輸出結(jié)果中都會包含，表示基因名字。

• log2FoldChange: 差異倍數(shù)的對數(shù)，一般的差異分析輸出結(jié)果中也會給出對數(shù)處理的值, 因此程序沒有提供這一步的計算操作。

• padj: 多重假設(shè)檢驗矯正過的差異顯著性P值；一般的差異分析輸出結(jié)果為原始值，程序提供一個參數(shù)對其求取負對數(shù)。

• significant: 可選列，標記哪些基因是上調(diào)、下調(diào)、無差異；若無此列或未在參數(shù)中指定此列，默認程序會根據(jù)padj列和log2FoldChange列根據(jù)給定的閾值自動計算差異基因，并作出不同顏色的標記。

• label?：可選列，一般用于在圖中標記出感興趣的基因的名字。非-行的字符串都會標記在圖上。

volcano = "id;log2FoldChange;padj;significant;label
E00007;4.28238;0;EHBIO_UP;A
E00008;-1.1036;0.476466843393901;Unchanged;-
E00009;-0.274368;1;Unchanged;-
E00010;4.62347;7.37606076333335e-103;EHBIO_UP;-
E00012;0.973987;0.482982440163204;Unchanged;-
E00017;-1.30205;0.000555693857439792;Baodian_UP;B
E00024;0.617636;2.78047837287061e-13;Unchanged;-
E00033;1.48669;2.56000581595275e-60;EHBIO_UP;-
E00034;-0.783716;0.00341521725291801;Unchanged;-
E00036;2.01592;6.03136656016401e-06;EHBIO_UP;C
E00040;-1.89657;4.73663890849056e-21;Baodian_UP;-
E00041;-0.268168;0.563429434558031;Unchanged;-
E00042;0.0861048;0.367700939634328;Unchanged;-
E00043;-1.19328;1.42673872027352e-153;Baodian_UP;-
E00044;-0.887981;2.43067804654905e-26;Unchanged;-
E00047;-0.610941;5.51696648645932e-57;Unchanged;-"

# 數(shù)據(jù)的讀取之前的R語言統(tǒng)計和繪圖系列都已解釋過，不再贅述
# 文末也有鏈接可直達之前的文章，新學者建議從頭開始
volcanoData <- read.table(text=volcano, sep=";", header=T, quote="", check.names=F)
head(volcanoData)

 ? ? ?id log2FoldChange ? ? ? ? ?padj significant label
1 E00007 ? ? ? 4.282380 ?0.000000e+00 ? ?EHBIO_UP ? ? A
2 E00008 ? ? ?-1.103600 ?4.764668e-01 ? Unchanged ? ? -
3 E00009 ? ? ?-0.274368 ?1.000000e+00 ? Unchanged ? ? -
4 E00010 ? ? ? 4.623470 7.376061e-103 ? ?EHBIO_UP ? ? -
5 E00012 ? ? ? 0.973987 ?4.829824e-01 ? Unchanged ? ? -
6 E00017 ? ? ?-1.302050 ?5.556939e-04 ?Baodian_UP ? ? B

繪制散點圖，只需要指定X軸和Y軸，再加上geom_point即可。

library(ggplot2)
p <- ggplot(volcanoData, aes(x=log2FoldChange, y=padj))
p <- p + geom_point()
# 前面是給p不斷添加圖層的過程
# 單輸入一個p是真正作圖
# 前面有人說，上面都輸完了，怎么沒出圖
# 就因為差了一個p
p

說好的火山圖的例子，但怎么也看不出噴發(fā)的態(tài)勢。

對數(shù)據(jù)坐下預處理，差異大的基因padj小，先對其求取負對數(shù)，所謂負負得正，差異大的基因就會處于圖的上方了。

# 從示例數(shù)據(jù)中看到，最小的padj值為0，求取負對數(shù)為正無窮。
# 實際上padj值小到一個點對我們來講就是個數(shù)
# 所以可以給所有小于1e-6的padj都讓其等于1e-6，再小也沒意義
#
volcanoData[volcanoData$padj<1e-6, "padj"] <- 1e-6
volcanoData$padj <- (-1)* log10(volcanoData$padj)

數(shù)據(jù)中基因的上調(diào)倍數(shù)遠高于下調(diào)倍數(shù)，使得出來的圖是偏的，這次畫圖時調(diào)整下X軸的區(qū)間使圖對稱。

p <- ggplot(volcanoData, ?aes(x=log2FoldChange, ?y=padj)) +
 ? ? geom_point() +
 ? ? xlim(-4.7, 4.7)
p

有點意思了，數(shù)據(jù)太少不明顯，下一步加上顏色看看。

p <- ggplot(volcanoData, ?aes(x=log2FoldChange, ?y=padj)) +
 ? ? geom_point(color=significant) +
 ? ? xlim(-4.7, 4.7)
p

利用現(xiàn)有的數(shù)據(jù)，基本上就是這個樣子了。雖然還不太像，原理都已經(jīng)都點到了。

Volcano plot | 別再問我這為什么是火山圖

盜取火山圖繪制一文中的圖來顯示個真正的火山圖吧。這樣一步步繪制很麻煩，去看一步法吧。

橫縱軸都為字符串的散點圖展示

輸入數(shù)據(jù)格式如下

這個數(shù)據(jù)是以前講到的FASTQC結(jié)果總結(jié)中的直觀的查看所有樣品測序堿基質(zhì)量和GC含量的散點圖的示例數(shù)據(jù)。

fastqc<-"ID;GC_quality;Base_quality
ehbio_1_1;PASS;PASS
ehbio_1_2;PASS;PASS
ehbio_2_1;WARN;PASS
ehbio_2_2;WARN;PASS
Other_1_1;FAIL;FAIL
Other_1_2;FAIL;FAIL"

fastqc_data <- read.table(text=fastqc, sep=";", header=T)
# 就不查看了

p <- ggplot(fastqc_data, aes(x=GC_quality, y=Base_quality)) + geom_point()
p

六個點少了只剩下了3個，重疊在一起了，而且也不知道哪個點代表什么樣品。這時需要把點抖動下，用到一個包ggbeeswarm，抖動圖的神器。

library(ggbeeswarm)
p <- ggplot(fastqc_data, aes(x=GC_quality, y=Base_quality)) + geom_quasirandom()
# 使用geom_text增加點的標記
# label表示標記哪一列的數(shù)值
# position_quasirandom獲取點偏移后的位置
# xjust調(diào)整對齊方式; hjust是水平的對齊方式，0為左，1為右，0.5居中，0-1之間可以取任意值。vjust是垂直對齊方式，0底對齊，1為頂對齊，0.5居中，0-1之間可以取任意值。
# check_overlap檢查名字在圖上是否重疊
p <- p + geom_text(aes(label=ID), position=position_quasirandom(),hjust=0, check_overlap=T)
p

一網(wǎng)打進散點圖繪制

假如有一個輸入數(shù)據(jù)如下所示(存儲于文件scatterplot.xls中)

Samp ? ?Gene1 ? ?Gene2 ? ?Color ? ?Size ? ?GC_quality ? ?Base_qualitya ? ?1 ? ?1 ? ?grp1 ? ?10 ? ?PASS ? ?PASSb ? ?2 ? ?2 ? ?grp1 ? ?10 ? ?PASS ? ?PASSc ? ?1 ? ?3 ? ?grp1 ? ?10 ? ?WARN ? ?PASSd ? ?3 ? ?1 ? ?grp2 ? ?15 ? ?WARN ? ?WARNe ? ?2 ? ?2 ? ?grp2 ? ?15 ? ?PASS ? ?WARNf ? ?3 ? ?3 ? ?grp3 ? ?5 ? ?PASS ? ?PASSg ? ?2 ? ?1 ? ?grp3 ? ?5 ? ?WARN ? ?PASS

想繪制樣品在這兩個Gene為軸的空間的分布，并標記樣品的屬性，只需要運行如下命令

# -f: 指定輸入文件，列數(shù)不限，順序不限; 第一行為列名字，第一列無特殊要求，必選
# -X: 指定哪一列為X軸信息，必選
# -Y: 指定哪一列為Y軸信息，必選
# -c: 指定用哪一列標記顏色，可選
# -s: 指定哪一列標記大小，一般為數(shù)字列，可選
# -S: 指定哪一列標記形狀，可選
# -L: 指定哪一列用來作為文本標記
# -w, -u: 指定圖的長寬
sp_scatterplot2.sh -f scatterplot.xls -X Gene1 -Y Gene2 -c Color -s Size -S GC_quality -L Samp -w 10 -u 10

如果橫縱軸為字符串，且有重復, 則需指定參數(shù)-J TRUE以錯開重疊的點，具體如下

# -O: 指定X軸變量的順序, 默認是字母順序
# 其它列或其它屬性的順序也可以用相應(yīng)的方式指示，具體看程序的幫助提示
# -c Gene1: 用特定基因的表達對點著色，單細胞分析圖中常用
# -J TRUE: 見上
# -Z FALSE：默認使用geom_text_repel添加點的標記，及其智能，不會出現(xiàn)標簽過多覆蓋的情況
# 但對jitterplot，會有些沖突，所以在`-J TRUE`且出來的圖中點的標簽不符合預期時，設(shè)定
# 次參數(shù)為FALSE，使用geom_text標記點。

sp_scatterplot2.sh -f scatterplot.xls -X GC_quality -Y Base_quality -O "'WARN', 'PASS'" -c Gene1 -w 10 -u 10 -J TRUE -L Samp -Z FALSE

只有想不到，沒有做不到，sp_scatterplot2.sh還可以完成更多你想做的散點圖，而且只需調(diào)參數(shù)，無需改代碼，簡單可重用。

 轉(zhuǎn)發(fā)朋友圈截圖給客服（shengxinbaodian）獲取，不用再重復了

Reference

• http://blog.genesino.com/2017/07/scatterPlot

• R語言 - 入門環(huán)境Rstudio

• R語言 - 熱圖繪制 (heatmap)

• R語言 - 基礎(chǔ)概念和矩陣操作

• R語言 - 熱圖簡化

• R語言 - 熱圖美化

• R語言 - 線圖繪制

• R語言 - 線圖一步法

• R語言 - 箱線圖（小提琴圖、抖動圖、區(qū)域散點圖）

• R語言 - 箱線圖一步法

• R語言 - 火山圖

高顏值免費在線繪圖（點擊圖片直達網(wǎng)站手冊）

往期精品(點擊圖片直達文字對應(yīng)教程)

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集