R语言学习 – 箱线图（小提琴图、抖动图、区域散点图）

来源： 生信宝典评论14,884

箱线图

箱线图是能同时反映数据统计量和整体分布，又很漂亮的展示图。在2014年的Nature Method上有2篇Correspondence论述了使用箱线图的好处和一个在线绘制箱线图的工具。就这样都可以发两篇Nature method，没天理，但也说明了箱线图的重要意义。

下面这张图展示了Bar plot、Box plot、Volin plot和Bean plot对数据分布的反应。从Bar plot上只能看到数据标准差或标准误不同；Box plot可以看到数据分布的集中性不同；Violin plot和Bean plot展示的是数据真正的分布，尤其是对Biomodal数据的展示。

Boxplot从下到上展示的是最小值，第一四分位数 (箱子的下边线)、中位数 (箱子中间的线)、第三四分位数 (箱子上边线)、最大值。

R语言学习 – 箱线图（小提琴图、抖动图、区域散点图）-图片1

http://www.nature.com/nmeth/journal/v11/n2/full/nmeth.2811.html

一步步解析箱线图绘制

假设有这么一个基因表达矩阵，第一列为基因名字，后面几列为样品名字，想绘制下样品中基因表达的整体分布。

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

读入数据并转换为ggplot2需要的长数据表格式 (经过前面几篇的练习，这应该都很熟了)

  1. profile_text <- read.table(text=profile, header=T, row.names=1, quote="",sep=";", check.names=F)
  2. # 在melt时保留位置信息
  3. # melt格式是ggplot2画图最喜欢的格式
  4. # 好好体会下这个格式，虽然多占用了不少空间，但是确实很方便
  5.  
  6. library(ggplot2)
  7. library(reshape2)
  8. data_m <- melt(profile_text)
  9. head(data_m)

 

  1.   variable value
  2. 1  2cell_1     4
  3. 2  2cell_1     6
  4. 3  2cell_1     8
  5. 4  2cell_1    10
  6. 5  2cell_1    12
  7. 6  2cell_1    14

 

像往常一样，就可以直接画图了。

  1. # variable和value为矩阵melt后的两列的名字，内部变量, variable代表了点线的属性，value代表对应的值。
  2. p <- ggplot(data_m, aes(x=variable, y=value),color=variable) + 
  3. geom_boxplot() + 
  4. theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
  5. theme(legend.position="none")
  6. p
  7. # 图会存储在当前目录的Rplots.pdf文件中，如果用Rstudio，可以不运行dev.off()
  8. dev.off()

 

箱线图出来了，看上去还可以，再加点色彩。

R语言学习 – 箱线图（小提琴图、抖动图、区域散点图）-图片2

  1. # variable和value为矩阵melt后的两列的名字，内部变量, variable代表了点线的属性，value代表对应的值。
  2. p <- ggplot(data_m, aes(x=variable, y=value),color=variable) + 
  3. geom_boxplot(aes(fill=factor(variable))) + 
  4. theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
  5. theme(legend.position="none")
  6. p
  7. # 图会存储在当前目录的Rplots.pdf文件中，如果用Rstudio，可以不运行dev.off()
  8. dev.off()

 

R语言学习 – 箱线图（小提琴图、抖动图、区域散点图）-图片3

再看看Violin plot

  1. # variable和value为矩阵melt后的两列的名字，内部变量, variable代表了点线的属性，value代表对应的值。
  2. p <- ggplot(data_m, aes(x=variable, y=value),color=variable) + 
  3. geom_violin(aes(fill=factor(variable))) + 
  4. theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
  5. theme(legend.position="none")
  6. p
  7. # 图会存储在当前目录的Rplots.pdf文件中，如果用Rstudio，可以不运行dev.off()
  8. dev.off()

 

R语言学习 – 箱线图（小提琴图、抖动图、区域散点图）-图片4

还有Jitter plot (这里使用的是ggbeeswarm包)

  1. library(ggbeeswarm)
  2. # 为了更好的效果，只保留其中一个样品的数据
  3. # grepl类似于Linux的grep命令，获取特定模式的字符串
  4.  
  5. data_m2 <- data_m[grepl("_3", data_m$variable),]
  6.  
  7. # variable和value为矩阵melt后的两列的名字，内部变量, variable代表了点线的属性，value代表对应的值。
  8. p <- ggplot(data_m2, aes(x=variable, y=value),color=variable) + 
  9. geom_quasirandom(aes(colour=factor(variable))) + 
  10. theme_bw() + theme(panel.grid.major = element_blank(),
  11. panel.grid.minor = element_blank(), legend.key=element_blank()) +
  12. theme(legend.position="none")
  13. # 也可以用geom_jitter(aes(colour=factor(variable)))代替geom_quasirandom(aes(colour=factor(variable)))
  14. # 但个人认为geom_quasirandom给出的结果更有特色
  15.  
  16. ggsave(p, filename="jitterplot.pdf", width=14, height=8, units=c("cm"))

 

R语言学习 – 箱线图（小提琴图、抖动图、区域散点图）-图片5

绘制单个基因 (A)的箱线图

为了更好的展示效果，下面的矩阵增加了样品数量和样品的分组信息。

  1. 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
  2. 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
  3. 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"
  4.  
  5. profile_text <- read.table(text=profile, header=T, row.names=1, quote="",sep=";", check.names=F)
  6.  
  7. data_m = data.frame(t(profile_text['A',]))
  8. data_m$sample = rownames(data_m)
  9. # 只挑选显示部分
  10. # grepl前面已经讲过用于匹配
  11. data_m[grepl('_[123]', data_m$sample),]

 

  1.            A   sample
  2. 2cell_1  4.0  2cell_1
  3. 2cell_2  6.0  2cell_2
  4. 2cell_3  7.0  2cell_3
  5. 4cell_1  3.2  4cell_1
  6. 4cell_2  5.2  4cell_2
  7. 4cell_3  5.6  4cell_3
  8. zygote_1 2.0 zygote_1
  9. zygote_2 4.0 zygote_2
  10. zygote_3 3.0 zygote_3

 

获得样品分组信息 (这个例子比较特殊，样品的分组信息就是样品名字下划线前面的部分)

  1. # 可以利用strsplit分割，取出其前面的字符串
  2. # R中复杂的输出结果多数以列表的形式体现，在之前的矩阵操作教程中
  3. # 提到过用str函数来查看复杂结果的结构，并从中获取信息
  4. group = unlist(lapply(strsplit(data_m$sample,"_"), function(x) x[1]))
  5. data_m$group = group
  6. data_m[grepl('_[123]', data_m$sample),]

 

  1.            A   sample  group
  2. 2cell_1  4.0  2cell_1  2cell
  3. 2cell_2  6.0  2cell_2  2cell
  4. 2cell_3  7.0  2cell_3  2cell
  5. 4cell_1  3.2  4cell_1  4cell
  6. 4cell_2  5.2  4cell_2  4cell
  7. 4cell_3  5.6  4cell_3  4cell
  8. zygote_1 2.0 zygote_1 zygote
  9. zygote_2 4.0 zygote_2 zygote
  10. zygote_3 3.0 zygote_3 zygote

 

如果没有这个规律，也可以提到类似于下面的文件，指定样品所属的组的信息。

  1. sampleGroup_text="Sample;Group
  2. zygote_1;zygote
  3. zygote_2;zygote
  4. zygote_3;zygote
  5. zygote_4;zygote
  6. zygote_5;zygote
  7. zygote_6;zygote
  8. 2cell_1;2cell
  9. 2cell_2;2cell
  10. 2cell_3;2cell
  11. 2cell_4;2cell
  12. 2cell_5;2cell
  13. 2cell_6;2cell
  14. 4cell_1;4cell
  15. 4cell_2;4cell
  16. 4cell_3;4cell
  17. 4cell_4;4cell
  18. 4cell_5;4cell
  19. 4cell_6;4cell"
  20.  
  21. #sampleGroup = read.table(text=sampleGroup_text,sep="\t",header=1,check.names=F,row.names=1)
  22.  
  23. #data_m <- merge(data_m, sampleGroup, by="row.names")
  24.  
  25. # 会获得相同的结果，脚本注释掉了以免重复执行引起问题。

 

矩阵准备好了，开始画图了 (小提琴图做例子，其它类似)

  1. # 调整下样品出现的顺序
  2. data_m$group <- factor(data_m$group, levels=c("zygote","2cell","4cell"))
  3. # group和A为矩阵中两列的名字，group代表了值的属性，A代表基因A对应的表达值。
  4. # 注意看修改了的地方
  5. p <- ggplot(data_m, aes(x=group, y=A),color=group) + 
  6. geom_violin(aes(fill=factor(group))) + 
  7. theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
  8. theme(legend.position="none")
  9. p
  10. # 图会存储在当前目录的Rplots.pdf文件中，如果用Rstudio，可以不运行dev.off()

 

R语言学习 – 箱线图（小提琴图、抖动图、区域散点图）-图片6

长矩阵绘制箱线图

常规矩阵绘制箱线图要求必须是个方正的矩阵输入，而有时想比较的几个组里面检测的值数目不同。比如有三个组，GrpA组检测了6个病人，GrpB组检测了10个病人，GrpC组是12个正常人的检测数据。这时就很难形成一个行位检测值，列为样品的矩阵，长表格模式就适合与这种情况。

  1. long_table <- "Grp;Value
  2. GrpA;10
  3. GrpA;11
  4. GrpA;12
  5. GrpB;5
  6. GrpB;4
  7. GrpB;3
  8. GrpB;2
  9. GrpC;2
  10. GrpC;3"
  11.  
  12. long_table <- read.table(text=long_table,sep="\t",header=1,check.names=F)
  13.  
  14. p <- ggplot(data_m, aes(x=Grp, y=Value),color=Grp) + 
  15. geom_violin(aes(fill=factor(Grp))) + 
  16. theme(axis.text.x=element_text(angle=50,hjust=0.5, vjust=0.5)) +
  17. theme(legend.position="none")
  18. p
  19. dev.off()

长表格形式自身就是常规矩阵melt后的格式，这种用来绘制箱线图就很简单了，就不做解释了。

发表评论

匿名网友

拖动滑块以完成验证
加载失败