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數據的展現。

 

Box plot從下到上展現的是最小值，第一四分位數 (箱子的下邊線)、中位數 (箱子中間的線)、第三四分位數 (箱子上邊線)、最大值，具體解讀看這裏 擴增子圖表解讀1箱線圖：Alpha多樣性。

一步步解析箱線圖繪製

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

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"

讀入數據並轉換爲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)
  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

像往常同樣，就能夠直接畫圖了。

# 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()

箱線圖出來了，看上去還能夠，再加點色彩

# 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()

再看看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()

還有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"))

繪製單個基因 (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),]

得到樣品分組信息 (這個例子比較特殊，樣品的分組信息就是樣品名字下劃線前面的部分)

# 能夠利用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),]

若是沒有這個規律，也能夠提到相似於下面的文件，指定樣品所屬的組的信息。

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")
# 會得到相同的結果，腳本註釋掉了以避免重複執行引發問題

矩陣準備好了，開始畫圖了 (小提琴圖作例子，其它相似)

# 調整下樣品出現的順序
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()

長矩陣繪製箱線圖

常規矩陣繪製箱線圖要求必須是個方正的矩陣輸入，而有時想比較的幾個組裏面檢測的值數目不一樣。好比有三個組，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(long_table, 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

長表格形式自身就是常規矩陣melt後的格式，這種用來繪製箱線圖就很簡單了，就不作解釋了。