生信寶典公眾號:http://mp.weixin.qq.com/s/c2nwATMwyU_9FgLprqmQ5A
散點圖
散點圖在生物信息分析中是應(yīng)用比較廣的一個圖衡怀,常見的差異基因火山圖扎瓶、功能富集分析泡泡圖、相關(guān)性分析散點圖码俩、抖動圖盖灸、PCA樣品分類圖等均驶。凡是想展示分布狀態(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是真正作圖
# 前面有人說般妙,上面都輸完了,怎么沒出圖
# 就因為差了一個pp
說好的火山圖的例子歪架,但怎么也看不出噴發(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-6volcanoData$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)都點到了稿饰。
盜取火山圖繪制一文中的圖來顯示個真正的火山圖吧锦秒。這樣一步步繪制很麻煩,去看一步法吧喉镰。
橫縱軸都為字符串的散點圖展示
輸入數(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_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為軸的空間的分布,并標記樣品的屬性撬讽,只需要運行如下命令
# -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ā)獲取牺六,不用再重復了
Reference
http://blog.genesino.com/2017/07/scatterPlot