简介
在查阅文献的过程中，看到了几幅非常不错的出版图，今天就跟着小编一起学习下，他们是怎么使用 R 绘制出来的。
今天主要介绍 第五幅图（E) —— 带置信区间的拟合曲线图及带误差项的散点图。这个图在科研绘图中较为常用，例如：绘制产品失效时间的散点图。
前四幅图的详细代码介绍可见：基于 R 语言的科研论文绘图技巧详解（4）、基于 R 语言的科研论文绘图技巧详解（3）基于 R 语言的科研论文绘图技巧详解（2）基于 R 语言的科研论文绘图技巧详解（1） 。最后一幅图会随后继续介绍，读者在学习过程中，可以将内部学到的知识点应用到自己的图形绘制中。推文已经将主要知识点进行罗列，更有利于读者学习和查阅。
那我们来看看，作者是怎么实现这个功能的吧，本文知识点较多，大家耐心学习，建议自己实践。对应代码、数据可在 GitHub - marco-meer/scifig_plot_examples_R: Scientific publication figure plotting examples with R[1] 中找到。
主要知识点
- 绘制散点图、丝带形状图；
- 绘制横向、纵向误差图；
- 学会小技巧：展示轴外部的图形。
绘图
加载包
首先加载一些需要使用到的包。
library(ggplot2) # Grammar of graphics
设置主题
接下来，为了方便起见，作者在绘图前设置好了主题，并将该函数命名为 my_theme。
这一部分在第一篇推文给出，代码将在文末中完整代码给出。
手动修改大部分面板，具体可以参考本篇文章[2]。或者观看我在 B 站发布的《R 语言可视化教程》，里面也有一些简单主题设置介绍。
导入数据
首先使用 read.csv() 导入三个基因相关数据，最后合并到 data_E 数据框中。前 6 行数据如下所示。
data_Ea = read.csv("./data_Ea.csv") data_Eb = read.csv("./data_Eb.csv") data_Ec = read.csv("./data_Ec.csv") data_E = data.frame("gene"=c(rep("gene a",nrow(data_Ea)), rep("gene b",nrow(data_Eb)), rep("gene c",nrow(data_Ec))), "t"=c(data_Ea$t,data_Eb$t,data_Ec$t), "C"=c(data_Ea$C,data_Eb$C,data_Ec$C), "pos_err"=c(data_Ea$err,data_Eb$pos_err,data_Ec$pos_err_C), "neg_err"=c(data_Ea$err,data_Eb$neg_err,data_Ec$neg_err_C), "pos_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$pos_err_t), "neg_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$neg_err_t) ) head(data_E) # gene t C pos_err neg_err pos_err_t neg_err_t # 1 gene a 12 0.27 0.03 0.03 NA NA # 2 gene a 24 0.19 0.02 0.02 NA NA # 3 gene a 36 0.17 0.01 0.01 NA NA # 4 gene a 48 0.13 0.04 0.04 NA NA # 5 gene a 60 0.09 0.03 0.03 NA NA # 6 gene a 72 0.10 0.01 0.01 NA NA
gene 是分类变量，t 表示时间将作为 x 轴。C 表示正常水平下的结果（小编对基因方面不是很懂，可能理解的不对）。pos_err 和 neg_err 是正向与反向误差。
自定义函数
首先，定义三个函数，以便稍后拟合散点图。这里的参数具体怎么得到的，作者没有介绍，小编不是很了解。或许是使用了某种参数估计的方法吧。
f1 = function(t) 0.2*exp(-t/48) f2 = function(t) 0.3*exp(-t/60) f3 = function(t) 0.4*exp(-t/72) t = seq(0,96,1) # 构造 f2 函数的数据 ribbon = data.frame( "f2" = 0.3*exp(-t/60), "t" = t )
绘制图形
定义图形形状分别为实心正方形、圆形、三角形 manual_pch =c(15,16,17) 。之后简单绘制下散点图geom_point()和 丝带形状的图 geom_ribbon()。
注意：这里的
factor()将 gene 变成离散形式。由于geom_ribbon()中使用的数据不同，所以要重新指定下新数据。
manual_pch =c(15,16,17) 定义图形形状 ggplot(data=data_E) + geom_point(aes(x=t,y=C,pch=factor(gene)),size=2) + geom_ribbon(data=ribbon, aes(x=t,ymin=0.85*f2,ymax=1.15*f2), fill="black",alpha=0.1)
在这个基础上使用 stat_function() 添加三个函数，两条虚线，一条实线。
stat_function(fun=f1, geom="line",linetype="dashed") + stat_function(fun=f2, geom="line") + stat_function(fun=f3, geom="line",linetype="dotted")
使用 geom_errorbar() 和 geom_errorbarh() 添加误差棒（纵向与横向）。注意内部参数不同：（ymin，ymax） 和 （xmin，xmax）。
geom_errorbar(aes(x=t,ymin=C-neg_err, ymax=C+pos_err), width=2) + geom_errorbarh(aes(y=C,xmin=t-neg_err_t,xmax=t+pos_err_t))
最后就是主题的设置了，包括运用了自定义的主题 my_theme()。
scale_shape_manual(values=manual_pch) + my_theme() + theme(legend.title=element_blank())+ theme(legend.position = c(0.9,0.95)) + scale_x_continuous(expand = c(0, 0), breaks = c(seq(0,96,12)), limits = c(0,96) ) + scale_y_continuous(expand = c(0, 0), breaks = c(seq(0,0.4,0.05)), limits = c(0,0.4) ) + theme( panel.grid.major = element_line("gray95", size = 0.1), # putting label closer to axis bc exponent makes it bigger axis.title.y = element_text(margin = margin(r = -9)) ) + xlab("time (h)") + ylab(expression(paste("concentration (mmol",~cm^-3,")")))
但是，最右边几个点并没有显示清楚，这里作者运用了一个小技巧，将这些点清晰的呈现出来。来看看这个代码的效果：
coord_cartesian(clip = "off") # to allow for plotting outside axes
完整代码
# Panel E ---- library(ggplot2) base_size = 12 my_theme <- function() { theme( aspect.ratio = 1, axis.line =element_line(colour = "black"), # shift axis text closer to axis bc ticks are facing inwards axis.text.x = element_text(size = base_size*0.8, color = "black", lineheight = 0.9, margin=unit(c(0.3,0.3,0.3,0.3), "cm")), axis.text.y = element_text(size = base_size*0.8, color = "black", lineheight = 0.9, margin=unit(c(0.3,0.3,0.3,0.3), "cm")), axis.ticks = element_line(color = "black", size = 0.2), axis.title.x = element_text(size = base_size, color = "black", margin = margin(t = -5)), # t (top), r (right), b (bottom), l (left) axis.title.y = element_text(size = base_size, color = "black", angle = 90, margin = margin(r = -5)), axis.ticks.length = unit(-0.3, "lines"), legend.background = element_rect(color = NA, fill = NA), legend.key = element_rect(color = "black", fill = "white"), legend.key.size = unit(0.5, "lines"), legend.key.height =NULL, legend.key.width = NULL, legend.text = element_text(size = 0.6*base_size, color = "black"), legend.title = element_text(size = 0.6*base_size, face = "bold", hjust = 0, color = "black"), legend.text.align = NULL, legend.title.align = NULL, legend.direction = "vertical", legend.box = NULL, panel.background = element_rect(fill = "white", color = NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), plot.title = element_text(size = base_size, color = "black"), ) } data_Ea = read.csv("./data_Ea.csv") data_Eb = read.csv("./data_Eb.csv") data_Ec = read.csv("./data_Ec.csv") data_E = data.frame("gene"=c(rep("gene a",nrow(data_Ea)), rep("gene b",nrow(data_Eb)), rep("gene c",nrow(data_Ec))), "t"=c(data_Ea$t,data_Eb$t,data_Ec$t), "C"=c(data_Ea$C,data_Eb$C,data_Ec$C), "pos_err"=c(data_Ea$err,data_Eb$pos_err,data_Ec$pos_err_C), "neg_err"=c(data_Ea$err,data_Eb$neg_err,data_Ec$neg_err_C), "pos_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$pos_err_t), "neg_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$neg_err_t) ) head(data_E) f1 = function(t) 0.2*exp(-t/48) f2 = function(t) 0.3*exp(-t/60) f3 = function(t) 0.4*exp(-t/72) t = seq(0,96,1) ribbon = data.frame( "f2" = 0.3*exp(-t/60), "t" = t ) head(ribbon) manual_pch =c(15,16,17) # available pch: type ?pch panel_E <- ggplot(data=data_E) + geom_point(aes(x=t,y=C,pch=factor(gene)),size=2) + geom_ribbon(data=ribbon, aes(x=t,ymin=0.85*f2,ymax=1.15*f2), fill="black",alpha=0.1) + stat_function(fun=f1, geom="line",linetype="dashed") + stat_function(fun=f2, geom="line") + stat_function(fun=f3, geom="line",linetype="dotted") + geom_errorbar(aes(x=t,ymin=C-neg_err, ymax=C+pos_err), width=2) + geom_errorbarh(aes(y=C,xmin=t-neg_err_t,xmax=t+pos_err_t))+ scale_shape_manual(values=manual_pch) + my_theme() + theme(legend.title=element_blank())+ theme(legend.position = c(0.9,0.95)) + scale_x_continuous(expand = c(0, 0), breaks = c(seq(0,96,12)), limits = c(0,96) ) + scale_y_continuous(expand = c(0, 0), breaks = c(seq(0,0.4,0.05)), limits = c(0,0.4) ) + theme( panel.grid.major = element_line("gray95", size = 0.1), # putting label closer to axis bc exponent makes it bigger axis.title.y = element_text(margin = margin(r = -9)) ) + xlab("time (h)") + ylab(expression(paste("concentration (mmol",~cm^-3,")"))) + coord_cartesian(clip = "off") # to allow for plotting outside axes panel_E
小编有话说
本文主要学到的知识点如下：
- 使用
geom_point()绘制散点图，geom_ribbon()绘制丝带形状图； - 使用
stat_function()添加函数曲线； - 使用
geom_errorbar()和geom_errorbarh添加误差棒（纵向与横向）； - 使用
coord_cartesian(clip = "off")允许展示外轴的图形。
