基于 R 语言的科研论文绘图技巧详解（5）
点击下方公众号，回复资料分享，收获惊喜
简介
在查阅文献的过程中，看到了几幅非常不错的出版图，今天就跟着小编一起学习下，他们是怎么使用 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")允许展示外轴的图形。
看完这篇文章，相信你的技能包又多了点东西。记得实操噢！看了不代表就会了～ 如果觉得内容有用的话，小编写的有心的话。给小编来杯咖啡吧！
参考资料
GitHub - marco-meer/scifig_plot_examples_R: Scientific publication figure plotting examples with R: https://github.com/marco-meer/scifig_plot_examples_R
[2]文章: https://ggplot2.tidyverse.org/reference/theme.html
推荐： 可以保存以下照片，在b站扫该二维码，或者b站搜索【庄闪闪】观看Rmarkdown系列的视频教程。Rmarkdown视频新增两节视频（写轮眼幻灯片制作）需要视频内的文档，可在公众号回复【rmarkdown】
R沟通｜Rmarkdown教程（4）
R沟通｜Rmarkdown教程（3）
R沟通｜Rmarkdown教程（2）