1. 程式人生 > >R學習箱線圖

R學習箱線圖



箱線圖

箱線圖是能同時反映資料統計量和整體分佈,又很漂亮的展示圖。在2014年的Nature Method上有2篇Correspondence論述了使用箱線圖的好處和一個線上繪製箱線圖的工具。就這樣都可以發兩篇Nature method,沒天理,但也說明了箱線圖的重要意義。

下面這張圖展示了Bar plot、Box plot、Volin plot和Bean plot對資料分佈的反應。從Bar plot上只能看到資料標準差或標準誤不同;Box plot可以看到資料分佈的集中性不同;Violin plot和Bean plot展示的是資料真正的分佈,尤其是對Biomodal資料的展示。

![]({{ site.img_url}}/splot/boxplot_nm.png)

一步步解析箱線圖繪製

假設有這麼一個基因表達矩陣,第一列為基因名字,後面幾列為樣品名字,想繪製下樣品中基因表達的整體分佈。

profile="Name;2cell_1;2cell_2;2cell_3;4cell_1;4cell_2;4cell_3;zygote_1;zygote_2;zygote_3
A;4;6;7;3.2;5.2;5.6;2;4;3
B;6;8;9;5.2;7.2;7.6;4;6;5
C;8;10;11;7.2;9.2;9.6;6;8;7
D;10;12;13;9.2;11.2;11.6;8;10;9
E;12;14;15;11.2;13.2;13.6;10;12;11
F;14;16;17;13.2;15.2;15.6;12;14;13
G;15;17;18;14.2;16.2;16.6;13;15;14
H;16;18;19;15.2;17.2;17.6;14;16;15
I;17;19;20;16.2;18.2;18.6;15;17;16
J;18;20;21;17.2;19.2;19.6;16;18;17
L;19;21;22;18.2;20.2;20.6;17;19;18
M;20;22;23;19.2;21.2;21.6;18;20;19
N;21;23;24;20.2;22.2;22.6;19;21;20
O;22;24;25;21.2;23.2;23.6;20;22;21"
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15

讀入資料並轉換為ggplot2需要的長資料表格式 (經過前面幾篇的練習,這應該都很熟了)

profile_text <- read.table(text=profile, header=T, row.names=1, quote="",sep=";", check.names=F)
# 在melt時保留位置資訊
# melt格式是ggplot2畫圖最喜歡的格式
# 好好體會下這個格式,雖然多佔用了不少空間,但是確實很方便

library(ggplot2)
library(reshape2)
data_m <- melt(profile_text)
head(data_m)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  variable value
1  2cell_1     4
2  2cell_1     6
3  2cell_1     8
4  2cell_1    10
5  2cell_1    12
6  2cell_1    14
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8

像往常一樣,就可以直接畫圖了。

# variable和value為矩陣melt後的兩列的名字,內部變數, variable代表了點線的屬性,value代表對應的值。
p <- ggplot(data_m, aes(x=variable, y=value),color=variable) + 
geom_boxplot() + 
theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
theme(legend.position="none")
p
# 圖會儲存在當前目錄的Rplots.pdf檔案中,如果用Rstudio,可以不執行dev.off()
dev.off()
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8

箱線圖出來了,看上去還可以,再加點色彩。

# variable和value為矩陣melt後的兩列的名字,內部變數, variable代表了點線的屬性,value代表對應的值。
p <- ggplot(data_m, aes(x=variable, y=value),color=variable) + 
geom_boxplot(aes(fill=factor(variable))) + 
theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
theme(legend.position="none")
p
# 圖會儲存在當前目錄的Rplots.pdf檔案中,如果用Rstudio,可以不執行dev.off()
dev.off()
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8

再看看Violin plot

# variable和value為矩陣melt後的兩列的名字,內部變數, variable代表了點線的屬性,value代表對應的值。
p <- ggplot(data_m, aes(x=variable, y=value),color=variable) + 
geom_violin(aes(fill=factor(variable))) + 
theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
theme(legend.position="none")
p
# 圖會儲存在當前目錄的Rplots.pdf檔案中,如果用Rstudio,可以不執行dev.off()
dev.off()
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8

還有Jitter plot (這裡使用的是ggbeeswarm包)

library(ggbeeswarm)
# 為了更好的效果,只保留其中一個樣品的資料
# grepl類似於Linux的grep命令,獲取特定模式的字串

data_m2 <- data_m[grepl("_3", data_m$variable),]

# variable和value為矩陣melt後的兩列的名字,內部變數, variable代表了點線的屬性,value代表對應的值。
p <- ggplot(data_m2, aes(x=variable, y=value),color=variable) + 
geom_quasirandom(aes(colour=factor(variable))) + 
theme_bw() + theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), legend.key=element_blank()) +
theme(legend.position="none")
# 也可以用geom_jitter(aes(colour=factor(variable)))代替geom_quasirandom(aes(colour=factor(variable)))
# 但個人認為geom_quasirandom給出的結果更有特色

ggsave(p, filename="jitterplot.pdf", width=14, height=8, units=c("cm"))
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16

繪製單個基因 (A)的箱線圖

為了更好的展示效果,下面的矩陣增加了樣品數量和樣品的分組資訊。

profile="Name;2cell_1;2cell_2;2cell_3;2cell_4;2cell_5;2cell_6;4cell_1;4cell_2;4cell_3;4cell_4;4cell_5;4cell_6;zygote_1;zygote_2;zygote_3;zygote_4;zygote_5;zygote_6
A;4;6;7;5;8;6;3.2;5.2;5.6;3.6;7.6;4.8;2;4;3;2;4;2.5
B;6;8;9;7;10;8;5.2;7.2;7.6;5.6;9.6;6.8;4;6;5;4;6;4.5"

profile_text <- read.table(text=profile, header=T, row.names=1, quote="",sep=";", check.names=F)

data_m = data.frame(t(profile_text['A',]))
data_m$sample = rownames(data_m)
# 只挑選顯示部分
# grepl前面已經講過用於匹配
data_m[grepl('_[123]', data_m$sample),]
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
           A   sample
2cell_1  4.0  2cell_1
2cell_2  6.0  2cell_2
2cell_3  7.0  2cell_3
4cell_1  3.2  4cell_1
4cell_2  5.2  4cell_2
4cell_3  5.6  4cell_3
zygote_1 2.0 zygote_1
zygote_2 4.0 zygote_2
zygote_3 3.0 zygote_3
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11

獲得樣品分組資訊 (這個例子比較特殊,樣品的分組資訊就是樣品名字下劃線前面的部分)

# 可以利用strsplit分割,取出其前面的字串
# R中複雜的輸出結果多數以列表的形式體現,在之前的矩陣操作教程中
# 提到過用str函式來檢視複雜結果的結構,並從中獲取資訊
group = unlist(lapply(strsplit(data_m$sample,"_"), function(x) x[1]))
data_m$group = group
data_m[grepl('_[123]', data_m$sample),]
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
           A   sample  group
2cell_1  4.0  2cell_1  2cell
2cell_2  6.0  2cell_2  2cell
2cell_3  7.0  2cell_3  2cell
4cell_1  3.2  4cell_1  4cell
4cell_2  5.2  4cell_2  4cell
4cell_3  5.6  4cell_3  4cell
zygote_1 2.0 zygote_1 zygote
zygote_2 4.0 zygote_2 zygote
zygote_3 3.0 zygote_3 zygote
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11

如果沒有這個規律,也可以提到類似於下面的檔案,指定樣品所屬的組的資訊。

sampleGroup_text="Sample;Group
zygote_1;zygote
zygote_2;zygote
zygote_3;zygote
zygote_4;zygote
zygote_5;zygote
zygote_6;zygote
2cell_1;2cell
2cell_2;2cell
2cell_3;2cell
2cell_4;2cell
2cell_5;2cell
2cell_6;2cell
4cell_1;4cell
4cell_2;4cell
4cell_3;4cell
4cell_4;4cell
4cell_5;4cell
4cell_6;4cell"

#sampleGroup = read.table(text=sampleGroup_text,sep="\t",header=1,check.names=F,row.names=1)

#data_m <- merge(data_m, sampleGroup, by="row.names")

# 會獲得相同的結果,指令碼註釋掉了以免重複執行引起問題。
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25

矩陣準備好了,開始畫圖了 (小提琴圖做例子,其它類似)

# 調整下樣品出現的順序
data_m$group <- factor(data_m$group, levels=c("zygote","2cell","4cell"))
# group和A為矩陣中兩列的名字,group代表了值的屬性,A代表基因A對應的表達值。
# 注意看修改了的地方
p <- ggplot(data_m, aes(x=group, y=A),color=group) + 
geom_violin(aes(fill=factor(group))) + 
theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
theme(legend.position="none")
p
# 圖會儲存在當前目錄的Rplots.pdf檔案中,如果用Rstudio,可以不執行dev.off()
dev.off()
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11

長矩陣繪製箱線圖

常規矩陣繪製箱線圖要求必須是個方正的矩陣輸入,而有時想比較的幾個組裡面檢測的值數目不同。比如有三個組,GrpA組檢測了6個病人,GrpB組檢測了10個病人,GrpC組是12個正常人的檢測資料。這時就很難形成一個行位檢測值,列為樣品的矩陣,長表格模式就適合與這種情況。

long_table <- "Grp;Value
GrpA;10
GrpA;11
GrpA;12
GrpB;5
GrpB;4
GrpB;3
GrpB;2
GrpC;2
GrpC;3"

long_table <- read.table(text=long_table,sep="\t",header=1,check.names=F)

p <- ggplot(data_m, aes(x=Grp, y=Value),color=Grp) + 
geom_violin(aes(fill=factor(Grp))) + 
theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
theme(legend.position="none")
p
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19

長表格形式自身就是常規矩陣melt後的格式,這種用來繪製箱線圖就很簡單了,就不舉例子了。

箱線圖 - 一步繪製

繪圖時通常會碰到兩個頭疼的問題:

  1. 有時需要繪製很多的圖,唯一的不同就是輸入檔案,其它都不需要修改。如果用R指令碼,需要反覆替換檔名,繁瑣又容易出錯。 (R也有命令列引數,不熟,有經驗的可以嘗試下)

  2. 每次繪圖都需要不斷的調整引數,時間久了不用,就忘記引數怎麼設定了;或者調整次數過多,有了很多版本,最後不知道用哪個了。

為了簡化繪圖、維持指令碼的一致,我用bash對繪圖命令做了一個封裝,通過配置修改命令列引數,生成相應的繪圖指令碼,然後再繪製。

首先把測試資料儲存到檔案中方便呼叫。資料矩陣儲存在boxplot.normal.datasampleGroupboxplot.melt.data檔案中 (TAB鍵分割,內容在文件最後。如果你手上有自己的資料,也可以拿來用)。

使用正常矩陣預設引數繪製箱線圖

# -f: 指定輸入的矩陣檔案,第一列為行名字,第一行為header
      列數不限,列名字不限;行數不限,行名字預設為文字
sp_boxplot.sh -f boxplot.normal.data
  • 1
  • 2
  • 3

箱線圖出來了,但有點小亂。

stat_smooth

# -f: 指定輸入的矩陣檔案,第一列為行名字,第一行為header
      列數不限,列名字不限;行數不限,行名字預設為文字
# -P: none, 去掉legend (uppercase P)
# -b: X-axis旋轉45度
# -V: TRUE 繪製小提琴圖
sp_boxplot.sh -f boxplot.normal.data -P none -b 45 -V TRUE
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6

stat_smooth

繪製單個基因的小提琴圖加抖動圖

# -q: 指定某一行的名字,此處為基因名,繪製基因A的表達圖譜
# -Q: 指定樣本分組,繪製基因A在不同樣品組的表達趨勢
# -F Group: sampleGroup中第二列的名字,指代分組資訊,根據需要修改
# -J TRUE: 繪製抖動圖 jitter plot
# -L: 設定X軸樣品組順序
# -c TRUE -C "'red', 'pink', 'blue'": 指定每個箱線圖的顏色
sp_boxplot.sh -f boxplot.normal.data -q A -Q sampleGroup -F Group -V TRUE -J TRUE -L "'zygote','2cell','4cell'" -c TRUE -C "'red', 'pink', 'blue'" -P none
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7

stat_smooth

使用melted矩陣預設引數繪箱線圖

# -f: 指定輸入檔案
# -m TRUE: 指定輸入的矩陣為melted format
# -d Expr:指定表達值所在的列
# -F Rep: 指定子類所在列,也就是legend 
# -a Group:指定X軸分組資訊
# -j TRUE: jitter plot
sp_boxplot.sh -f boxplot.melt.data -m TRUE -d Expr -F Rep -a Group  -j TRUE
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7

stat_smooth

# 如果沒有子類,則-a和-F指定為同一值
# -R TRUE: 旋轉boxplot
sp_boxplot.sh -f boxplot.melt.data -m TRUE -d Expr -a Group -F Group -J TRUE -R TRUE
  • 1
  • 2
  • 3

stat_smooth

引數中最需要注意的是引號的使用:

  • 外層引號與內層引號不能相同
  • 凡引數值中包括了空格括號逗號等都用引號括起來作為一個整體。

為了推廣,也為了激起大家的熱情,如果想要sp_boxplot.sh指令碼的,還需要勞煩大家動動手,轉發此文章到朋友圈,並留言索取。

也希望大家能一起開發,完善功能。

#boxplot.normal.data
Name    2cell_1 2cell_2 2cell_3 2cell_4 2cell_5 2cell_6 4cell_1 4cell_2 4cell_3 4cell_4 4cell_5 4cell_6 zygote_1    zygote_2    zygote_3    zygote_4    zygote_5    zygote_6
A   4   6   7   5   8   6   3.2 5.2 5.6 3.6 7.6 4.8 2   4   3   2   4   2.5
B   6   8   9   7   10  8   5.2 7.2 7.6 5.6 9.6 6.8 4   6   5   4   6   4.5
C   8   10  11  9   12  10  7.2 9.2 9.6 7.6 11.6    8.8 6   8   7   6   8   6.5
D   10  12  13  11  14  12  9.2 11.2    11.6    9.6 13.6    10.8    8   10  9   8   10  8.5
E   12  14  15  13  16  14  11.2    13.2    13.6    11.6    15.6    12.8    10  12  11  10  12  10.5
F   14  16  17  15  18  16  13.2    15.2    15.6    13.6    17.6    14.8    12  14  13  12  14  12.5
G   15  17  18  16  19  17  14.2    16.2    16.6    14.6    18.6    15.8    13  15  14  13  15  13.5
H   16  18  19  17  20  18  15.2    17.2    17.6    15.6    19.6    16.8    14  16  15  14  16  14.5
I   17  19  20  18  21  19  16.2    18.2    18.6    16.6    20.6    17.8    15  17  16  15  17  15.5
J   18  20  21  19  22  20  17.2    19.2    19.6    17.6    21.6    18.8    16  18  17  16  18  16.5
L   19  21  22  20  23  21  18.2    20.2    20.6    18.6    22.6    19.8    17  19  18  17  19  17.5
M   20  22  23  21  24  22  19.2    21.2    21.6    19.6    23.6    20.8    18  20  19  18  20  18.5
N   21  23  24  22  25  23  20.2    22.2    22.6    20.6    24.6    21.8    19  21  20  19  21  19.5
O   22  24  25  23  26  24  21.2    23.2    23.6    21.6    25.6    22.8    20  22  21  20  22  20.5
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
#boxplot.melt.data

Gene    Sample  Group   Expr    Rep
A   zygote_1    zygote  2   1
A   zygote_2    zygote  4   2
A   zygote_3    zygote  3   3
A   zygote_4    zygote  2   4
A   zygote_5    zygote  4   5
A   zygote_6    zygote  2.5 6
A   2cell_1 2cell   4   1
A   2cell_2 2cell   6   2
A   2cell_3 2cell   7   3
A   2cell_4 2cell   5   4
A   2cell_5 2cell   8   5
A   2cell_6 2cell   6   6
A   4cell_1 4cell   3.2 1
A   4cell_2 4cell   5.2 2
A   4cell_3 4cell   5.6 3
A   4cell_4 4cell   3.6 4
A   4cell_5 4cell   7.6 5
A   4cell_6 4cell   4.8 6
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
#sampleGroup
Sample  Group
zygote_1    zygote
zygote_2    zygote
zygote_3    zygote
zygote_4    zygote
zygote_5    zygote
zygote_6    zygote
2cell_1 2cell
2cell_2 2cell
2cell_3 2cell
2cell_4 2cell
2cell_5 2cell
2cell_6 2cell
4cell_1 4cell
4cell_2 4cell
4cell_3 4cell
4cell_4 4cell
4cell_5 4cell
4cell_6 4cell