散點(diǎn)圖

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

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

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

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

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

• id: 不是必須的，但一般的軟件輸出結(jié)果中都會包含，表示基因名字。
• log2FoldChange: 差異倍數(shù)的對數(shù)，一般的差異分析輸出結(jié)果中也會給出對數(shù)處理的值, 因此程序沒有提供這一步的計算操作。
• padj: 多重假設(shè)檢驗(yàn)矯正過的差異顯著性P值；一般的差異分析輸出結(jié)果為原始值，程序提供一個參數(shù)對其求取負(fù)對數(shù)。
• significant: 可選列，標(biāo)記哪些基因是上調(diào)、下調(diào)、無差異；若無此列或未在參數(shù)中指定此列，默認(rèn)程序會根據(jù)padj列和log2FoldChange列根據(jù)給定的閾值自動計算差異基因，并作出不同顏色的標(biāo)記。
• label: 可選列，一般用于在圖中標(biāo)記出感興趣的基因的名字。非-行的字符串都會標(biāo)記在圖上。

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)計和繪圖系列都已解釋過，不再贅述 # 文末也有鏈接可直達(dá)之前的文章，新學(xué)者建議從頭開始 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

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

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

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

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

# 從示例數(shù)據(jù)中看到，最小的padj值為0，求取負(fù)對數(shù)為正無窮。 # 實(shí)際上padj值小到一個點(diǎn)對我們來講就是個數(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ù)遠(yuǎn)高于下調(diào)倍數(shù)，使得出來的圖是偏的，這次畫圖時調(diào)整下X軸的區(qū)間使圖對稱。

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

有點(diǎn)意思了，數(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)都點(diǎn)到了。

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

橫縱軸都為字符串的散點(diǎn)圖展示

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

這個數(shù)據(jù)是前面講到的FASTQC結(jié)果總結(jié)中的直觀的查看所有樣品測序堿基質(zhì)量和GC含量的散點(diǎn)圖的示例數(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

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

library(ggbeeswarm) p <- ggplot(fastqc_data, aes(x=GC_quality, y=Base_quality)) + geom_quasirandom() # 使用geom_text增加點(diǎn)的標(biāo)記 # label表示標(biāo)記哪一列的數(shù)值 # position_quasirandom獲取點(diǎn)偏移后的位置 # 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)打進(jìn)散點(diǎn)圖繪制

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

Samp Gene1 Gene2 Color Size GC_quality Base_quality a 1 1 grp1 10 PASS PASS b 2 2 grp1 10 PASS PASS c 1 3 grp1 10 WARN PASS d 3 1 grp2 15 WARN WARN e 2 2 grp2 15 PASS WARN f 3 3 grp3 5 PASS PASS g 2 1 grp3 5 WARN PASS

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

# -f: 指定輸入文件，列數(shù)不限，順序不限; 第一行為列名字，第一列無特殊要求，必選 # -X: 指定哪一列為X軸信息，必選 # -Y: 指定哪一列為Y軸信息，必選 # -c: 指定用哪一列標(biāo)記顏色，可選 # -s: 指定哪一列標(biāo)記大小，一般為數(shù)字列，可選 # -S: 指定哪一列標(biāo)記形狀，可選 # -L: 指定哪一列用來作為文本標(biāo)記 # -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

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

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

轉(zhuǎn)發(fā)獲取，不用再重復(fù)了

Reference

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

生信寶典公眾號, 幾千人一起學(xué)生信

https://mp.weixin.qq.com/s/c2nwATMwyU_9FgLprqmQ5A

總結(jié)

以上是生活随笔為你收集整理的R 学习 - 散点图的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。