做科研的Graphpad Prism是避不开的一道关，那么是否可以用R来实现其绘图风格呢？在微信公众号例竟然发现一篇好文，《一行代码绘制Graphpad Prism风格学术表》，下面就来一起学习下吧，接下来是对这篇推文的记录和整理，原文链接在文末：

乍一看，该R包的图标就有浓浓的Graphpad风。

核心函数：

• theme_prism()：绘图主题

• scale_colour/fill/shape_prism()：颜色、填充、形状

• scale_x/y_discrete/continuous(guide = “prism_bracket”)：刻度设置 或者guides(y/x=“prism_bracket”)

• add_pvalues()：自动添加p值

安装

library(pacman)
p_load(ggprism)
p_load(tidyverse)

#开发者版本
#remotes::install_github("csdaw/ggprism")

快速入门

tg <- ToothGrowth
tg$dose <- as.factor(tg$dose)
base <- ggplot(tg,aes(dose,len))+
  geom_violin(aes(color=dose,fill=dose),trim=F)+
  geom_boxplot(aes(fill=dose),width=0.2,color="black")+
  scale_y_continuous(limits = c(-5,40))

p_vals <- tibble::tribble(
  ~group1, ~group2, ~p.adj,   ~y.position,
  "0.5",   "1",     8.80e-14, 35,
  "0.5",   "2",     1.27e-7,  39
)
p1 <- base+
  scale_color_prism("floral")+
  scale_fill_prism("floral")+
  guides(y="prism_offset_minor")+
  theme_prism(base_size=16)+
  theme(legend.position = "none")+
  add_pvalue(p_vals,label = "p = {p.adj}", tip.length = 0,label.size = 4)

#快速拼图
library(patchwork)
base+p1

image.png

p值添加

考虑到科研中大家最常关注的就是p值，下面把关注点着重放在pvalue上面：

模板一

p_load(rstatix)
df_p_val <- rstatix::t_test(tg,len~supp) %>% 
  add_x_position()
p <- ggplot(tg, aes(factor(supp),len))+
  stat_summary(geom = "col",fun=mean)+
  stat_summary(geom = "errorbar",
               fun=mean,
               fun.min = function(x) mean(x)-sd(x),
               fun.max = function(x) mean(x)+sd(x),
               width=0.3)+
  theme_prism()+
  coord_cartesian(ylim=c(0,35))+
  scale_y_continuous(breaks = seq(0,35,5),expand = c(0,0))
p+add_pvalue(df_p_val,y.position =30)

image.png

模板二

df_p_val <- rstatix::t_test(tg,len~dose,ref.group = "0.5") %>% 
  add_xy_position()
p <- ggplot(tg, aes(factor(dose),len))+
  stat_summary(geom = "col",fun=mean)+
  stat_summary(geom = "errorbar",
               fun=mean,
               fun.min = function(x) mean(x)-sd(x),
               fun.max = function(x) mean(x)+sd(x),
               width=0.3)+
  theme_prism()+
  coord_cartesian(ylim=c(0,40))+
  scale_y_continuous(breaks = seq(0,40,5),expand = c(0,0))
p1 <- p+add_pvalue(df_p_val,label = "p.adj.signif")
p2 <- p+add_pvalue(df_p_val,label = "p.adj.signif",remove.bracket = T)
p1+p2

image.png

进阶版

示例数据采用的是Prism8中XY的剂量反应关系。

rm(list = ls())
p_load(tidyverse)
p_load(ggprism)
p_load(ggnewscale)

# 输入数据准备

df <- data.frame(
  agonist = c(1e-10, 1e-8, 3e-8, 1e-7, 3e-7, 1e-6, 3e-6, 1e-5, 3e-5, 1e-4, 3e-4),
  ctr1 = c(0, 11, 125, 190, 258, 322, 354, 348, NA, 412, NA),
  ctr2 = c(3, 33, 141, 218, 289, 353, 359, 298, NA, 378, NA),
  ctr3 = c(2, 25, 160, 196, 345, 328, 369, 372, NA, 399, NA),
  trt1 = c(3, NA, 11, 52, 80, 171, 289, 272, 359, 352, 389),
  trt2 = c(5, NA, 25, 55, 77, 195, 230, 333, 306, 320, 338), 
  trt3 = c(4, NA, 28, 61, 44, 246, 243, 310, 297, 365, NA)
) %>% 
  mutate(log.agonist = log10(agonist)) %>% 
  pivot_longer(
    c(-agonist, -log.agonist), 
    names_pattern = "(.{3})([0-9])", 
    names_to = c("treatment", "rep"),
    values_to = "response"
  ) %>% 
  filter(!is.na(response))
head(df)

## # A tibble: 6 x 5
##        agonist log.agonist treatment rep   response
##          <dbl>       <dbl> <chr>     <chr>    <dbl>
## 1 0.0000000001         -10 ctr       1            0
## 2 0.0000000001         -10 ctr       2            3
## 3 0.0000000001         -10 ctr       3            2
## 4 0.0000000001         -10 trt       1            3
## 5 0.0000000001         -10 trt       2            5
## 6 0.0000000001         -10 trt       3            4

p <- ggplot(df,aes(x=log.agonist, y=response))
dose_resp <- y ~ min + ((max - min) / (1 + exp(hill_coefficient * (ec50 - x))))
p <- p+geom_smooth(aes(color=treatment),
                   method = "nls",
                   formula = dose_resp,
                   se=F,
                   method.args=list(start=list(min=1.67,max=397,ec50=-7,hill_coefficient=1)))
p

image.png

# 改变曲线颜色

p <- p+scale_color_manual(
  labels = c("No inhibitor", "Inhibitor"),
  values = c("#00167B", "#9FA3FE")
)
p

image.png

# 调整线条颜色
p <- p+ggnewscale::new_scale_color()+
  geom_point(aes(color=treatment,shape=treatment),size=3)+
  scale_color_prism(
    palette = "winter_bright",
    labels=c("No inhibitor",
             "Inhibitor"))+
  scale_shape_prism(
    labels=c("No inhibitor",
             "Inhibitor"))
p

image.png

# 调整主题
p <- p+theme_prism(palette = "winter_bright",
                   base_size = 16)
p

image.png

# 调整y轴
p <- p+scale_y_continuous(
  limits=c(-100,500),
  breaks = seq(-100,500,100),
  guide = "prism_offset"
)
p

image.png

#调整x轴
p <- p + scale_x_continuous(
  limits = c(-10, -3), 
  breaks = -10:-3,
  guide = "prism_offset_minor",
  minor_breaks = log10(rep(1:9, 7)*(10^rep(-10:-4, each = 9))),
  labels = function(lab) {
    do.call(
      expression,
      lapply(paste(lab), function(x) bquote(bold("10"^.(x))))
    )
  }
)
p

image.png

# 调整x轴和y轴标题，删除legend
p <- p+theme(
  axis.title.y = element_blank(),
  axis.title.x = element_blank(),
  legend.position = c(0.05,0.95),
  legend.justification = c(0.05,0.95)
)+
  labs(x="[Agonist], M")
p

image.png

完整代码(remember the order of layers is important)

dose_resp <- y ~ min + ((max - min) / (1 + exp(hill_coefficient * (ec50 - x))))

ggplot(df, aes(x = log.agonist, y = response)) + 
  geom_smooth(
    aes(colour = treatment),
    method = "nls", formula = dose_resp, se = FALSE,
    method.args = list(start = list(min = 1.67, max = 397, ec50 = -7, hill_coefficient = 1))
  ) + 
  scale_colour_manual(labels = c("No inhibitor", "Inhibitor"),
                      values = c("#00167B", "#9FA3FE")) + 
  ggnewscale::new_scale_colour() +
  geom_point(aes(colour = treatment, shape = treatment), size = 3) + 
  scale_colour_prism(palette = "winter_bright", 
                     labels = c("No inhibitor", "Inhibitor")) + 
  scale_shape_prism(labels = c("No inhibitor", "Inhibitor")) + 
  theme_prism(palette = "winter_bright", base_size = 16) + 
  scale_y_continuous(limits = c(-100, 500), 
                     breaks = seq(-100, 500, 100),
                     guide = "prism_offset") + 
  scale_x_continuous(
    limits = c(-10, -3), 
    breaks = -10:-3,
    guide = "prism_offset_minor",
    minor_breaks = log10(rep(1:9, 7)*(10^rep(-10:-4, each = 9))),
    labels = function(lab) {
      do.call(
        expression,
        lapply(paste(lab), function(x) bquote(bold("10"^.(x))))
      )
    }
  ) + 
  theme(axis.title.y = element_blank(),
        legend.title = element_blank(),
        legend.position = c(0.05, 0.95),
        legend.justification = c(0.05, 0.95)) + 
  labs(x = "[Agonist], M")

参考链接：

• ggprism

• 一行代码绘制Graphpad Prism风格学术图表