基于 R 语言的科研论文绘图技巧详解(3)

2022-05-24 15:45:17 浏览数 (1)

简介

在查阅文献的过程中,看到了几幅非常不错的出版图,今天就跟着小编一起学习下,他们是怎么使用 R 绘制出来的。

今天主要介绍 第三幅图(C)——两组数据的带抖动点箱线图,这种图形在小编的研究方向中,经常会用来比较不同参数估计方法(极大似然,贝叶斯方法等)的估计性能好坏。这幅图比较常规,难点在于前期的数据准备。

前两幅图的详细代码介绍可见:基于 R 语言的科研论文绘图技巧详解(2),基于 R 语言的科研论文绘图技巧详解(1)。后面几幅图会一一介绍,读者在学习过程中,可以将内部学到的知识点应用到自己的图形绘制中。推文已经将主要知识点进行罗列,更有利于读者学习和查阅。

那我们来看看,作者是怎么实现这个功能的吧,本文知识点较多,大家耐心学习,建议自己实践。对应代码、数据可在 GitHub - marco-meer/scifig_plot_examples_R: Scientific publication figure plotting examples with R[1] 中找到。

主要知识点

  • 学会转化数据为图形所需的数据格式;
  • 学会绘制三变量的箱线图;
  • 学会绘制带抖动的散点图并修改透明度。

绘图

加载包

首先加载一些需要使用到的包。

代码语言:javascript复制
library(ggplot2) # Grammar of graphics

设置主题

接下来,为了方便起见,作者在绘图前设置好了主题,并将该函数命名为 my_theme

这一部分在第一篇推文[]给出,代码将在文末中完整代码给出。

手动修改大部分面板,具体可以参考本篇文章[2]。或者观看我在 B 站发布的《R 语言可视化教程》,里面也有一些简单主题设置介绍。

导入数据

首先使用 read.csv() 导入数据,其中一个数据前几行如下所示。

代码语言:javascript复制
data_Cwt_E8.5  = read.csv( "./data_Cwt_E8.5.csv")
data_Cwt_E9.5  = read.csv( "./data_Cwt_E9.5.csv")
data_Cwt_E10.5 = read.csv("./data_Cwt_E10.5.csv")
data_Cwt_E11.5 = read.csv("./data_Cwt_E11.5.csv")
data_Cmu_E8.5  = read.csv( "./data_Cmu_E8.5.csv")
data_Cmu_E9.5  = read.csv( "./data_Cmu_E9.5.csv")
data_Cmu_E10.5 = read.csv("./data_Cmu_E10.5.csv")
data_Cmu_E11.5 = read.csv("./data_Cmu_E11.5.csv")
)
head(data_Cwt_E8.5)
#   X42.0231534315441
# 1          48.60862
# 2          59.44292
# 3          44.47093
# 4          45.87919
# 5          37.43912
# 6          54.70787

之后,使用data.frame()转化数据成 ggplot 所需要的格式。这里作者使用基础包中的 rep() 一列列的构造数据。当然如果你会使用 tidyverse,可以换种方式整理数据。具体可以参考:《R语言教程》[3]R 数据科学[4],或者公众号相应的数据科学系列文章。

代码语言:javascript复制
data_C = data.frame(
  "type" = c(
    rep(
      "wildtype",
      nrow(data_Cwt_E8.5)  
        nrow(data_Cwt_E9.5)  
        nrow(data_Cwt_E10.5)  
        nrow(data_Cwt_E11.
    ),
    rep(
      "mutant",
      nrow(data_Cmu_E8.5)  
        nrow(data_Cmu_E9.5)  
        nrow(data_Cmu_E10.5)  
        nrow(data_Cmu_E11.5)
    )
  ),
  "dev_stage" = c(
    rep("E8.5", nrow(data_Cwt_E8.5)),
    rep("E9.5", nrow(data_Cwt_E9.5)),
    rep("E10.5", nrow(data_Cwt_E10.5)),
    rep("E11.5", nrow(data_Cwt_E11.5)),
    rep("E8.5", nrow(data_Cmu_E8.5)),
    rep("E9.5", nrow(data_Cmu_E9.5)),
    rep("E10.5", nrow(data_Cmu_E10.5)),
    rep("E11.5", nrow(data_Cmu_E11.5))
  ),
  "trachea_length" = c(
    data_Cwt_E8.5[, 1] ,
    data_Cwt_E9.5[, 1] ,
    data_Cwt_E10.5[, 1] ,
    data_Cwt_E11.5[, 1] ,
    data_Cmu_E8.5[, 1] ,
    data_Cmu_E9.5[, 1] ,
    data_Cmu_E10.5[, 1] ,
    data_Cmu_E11.5[, 1]
  )
)

head(data_C)
#      type dev_stage trachea_length
# 1 wildtype      E8.5       48.60862
# 2 wildtype      E8.5       59.44292
# 3 wildtype      E8.5       44.47093
# 4 wildtype      E8.5       45.87919
# 5 wildtype      E8.5       37.43912
# 6 wildtype      E8.5       54.70787

这里得到的数据,一共有三列,不同的数据集的数据值在 trachea_length 中,typedev_stage 为离散数据。

绘图步骤详解

这幅图的绘图代码比较传统,但是还是有些细节需要和大家分享下。

基础版本

使用 geom_boxplot() 绘制箱线图,geom_point() 加入散点图,注意这里使用了 position = position_jitterdodge() 使得图形相对分散开,使用 alpha 修改图形透明度。

代码语言:javascript复制
  ggplot(data=data_C,
         aes(x=dev_stage,
             y=trachea_length,
             fill=type)) 
  geom_boxplot()   
  geom_point(position=position_jitterdodge(), # jitter for h-dist, dodge for grouped dists
             pch=21, #圆形
             alpha=0.4) 

可以看到,这幅图基本已经完成,为了提升美观,和图片的可读性,我们又做了后续的处理。

微调主题

首先,修改 x 轴标签 (scale_x_discrete),填充颜色配色(scale_fill_manual)。

代码语言:javascript复制
  scale_x_discrete(limits=c("E8.5","E9.5","E10.5","E11.5"))  
  scale_fill_manual(values=c('#f2a340','#998fc2')) 

之后添加主题,使用先前设定好的主题函数 my_theme()与其他细节调整。

my_theme() 的代码,在文末完整版代码中会给出。

代码语言:javascript复制
my_theme()  
  theme(legend.position = c(0.18,0.95),
        legend.title = element_blank(),
        legend.key = element_blank()
  )  
  scale_y_continuous(expand = c(0, 0),
                     breaks=c(seq(20,120,by=20)),
                     limits = c(20,120)
  ) 
  xlab("developmental stage (days)")   
  ylab(expression(paste("tracheal length (",~mu*m,")")))

完整代码

代码语言:javascript复制
# Panel C ----
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_Cwt_E8.5  = read.csv( "./data_Cwt_E8.5.csv")
data_Cwt_E9.5  = read.csv( "./data_Cwt_E9.5.csv")
data_Cwt_E10.5 = read.csv("./data_Cwt_E10.5.csv")
data_Cwt_E11.5 = read.csv("./data_Cwt_E11.5.csv")
data_Cmu_E8.5  = read.csv( "./data_Cmu_E8.5.csv")
data_Cmu_E9.5  = read.csv( "./data_Cmu_E9.5.csv")
data_Cmu_E10.5 = read.csv("./data_Cmu_E10.5.csv")
data_Cmu_E11.5 = read.csv("./data_Cmu_E11.5.csv")

head(data_Cwt_E8.5)
# format data to ggplot's liking
data_C = data.frame(
  "type" = c(
    rep(
      "wildtype",
      nrow(data_Cwt_E8.5)  
        nrow(data_Cwt_E9.5)  
        nrow(data_Cwt_E10.5)  
        nrow(data_Cwt_E11.5)
    ),
    rep(
      "mutant",
      nrow(data_Cmu_E8.5)  
        nrow(data_Cmu_E9.5)  
        nrow(data_Cmu_E10.5)  
        nrow(data_Cmu_E11.5)
    )
  ),
  "dev_stage" = c(
    rep("E8.5", nrow(data_Cwt_E8.5)),
    rep("E9.5", nrow(data_Cwt_E9.5)),
    rep("E10.5", nrow(data_Cwt_E10.5)),
    rep("E11.5", nrow(data_Cwt_E11.5)),
    rep("E8.5", nrow(data_Cmu_E8.5)),
    rep("E9.5", nrow(data_Cmu_E9.5)),
    rep("E10.5", nrow(data_Cmu_E10.5)),
    rep("E11.5", nrow(data_Cmu_E11.5))
  ),
  "trachea_length" = c(
    data_Cwt_E8.5[, 1] ,
    data_Cwt_E9.5[, 1] ,
    data_Cwt_E10.5[, 1] ,
    data_Cwt_E11.5[, 1] ,
    data_Cmu_E8.5[, 1] ,
    data_Cmu_E9.5[, 1] ,
    data_Cmu_E10.5[, 1] ,
    data_Cmu_E11.5[, 1]
  )
)

head(data_C)

panel_C <- 
  ggplot(data=data_C,
         aes(x=dev_stage,
             y=trachea_length,
             fill=type)) 
  geom_boxplot()   
  geom_point(position=position_jitterdodge(), # jitter for h-dist, dodge for grouped dists
             pch=21, #圆形
             alpha=0.4)    # transparency
  scale_x_discrete(limits=c("E8.5","E9.5","E10.5","E11.5"))  
  scale_fill_manual(values=c('#f2a340','#998fc2'))   # custom colors in hex code
  my_theme()  
  theme(legend.position = c(0.18,0.95),
        legend.title = element_blank(),
        legend.key = element_blank()
  )  
  scale_y_continuous(expand = c(0, 0),
                     breaks=c(seq(20,120,by=20)),
                     limits = c(20,120)
  ) 
  xlab("developmental stage (days)")   
  ylab(expression(paste("tracheal length (",~mu*m,")")))

panel_C

小编有话说

本文主要学到的知识点如下:

  1. 使用基础包的 data.frame()rep() 整理和转化数据;
  2. 使用 geom_boxplot() 绘制箱线图并添加第三个变量;
  3. 使用 position = position_jitterdodge() 将散点分散展示。

参考资料

[1]

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

[3]

《R语言教程》: https://www.math.pku.edu.cn/teachers/lidf/docs/Rbook/html/_Rbook/ggplotvis.html

[4]

R 数据科学: https://r4ds.had.co.nz/

0 人点赞