做科研的Graphpad Prism是避不開的一道關(guān)炎码，那么是否可以用R來實現(xiàn)其繪圖風(fēng)格呢？在微信公眾號例竟然發(fā)現(xiàn)一篇好文秋泳，《一行代碼繪制Graphpad Prism風(fēng)格學(xué)術(shù)表》潦闲，下面就來一起學(xué)習(xí)下吧，接下來是對這篇推文的記錄和整理迫皱，原文鏈接在文末：
乍一看歉闰，該R包的圖標(biāo)就有濃濃的Graphpad風(fēng)辖众。
核心函數(shù)：
theme_prism()：繪圖主題scale_colour/fill/shape_prism()：顏色、填充和敬、形狀scale_x/y_discrete/continuous(guide = “prism_bracket”)：刻度設(shè)置 或者guides(y/x=“prism_bracket”)add_pvalues()：自動添加p值
安裝
library(pacman)
p_load(ggprism)
p_load(tidyverse)
#開發(fā)者版本
#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值添加
考慮到科研中大家最常關(guān)注的就是p值凹炸，下面把關(guān)注點著重放在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
進階版
示例數(shù)據(jù)采用的是Prism8中XY的劑量反應(yīng)關(guān)系。
rm(list = ls())
p_load(tidyverse)
p_load(ggprism)
p_load(ggnewscale)
# 輸入數(shù)據(jù)準(zhǔn)備
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
# 調(diào)整線條顏色
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
# 調(diào)整主題
p <- p+theme_prism(palette = "winter_bright",
base_size = 16)
p
image.png
# 調(diào)整y軸
p <- p+scale_y_continuous(
limits=c(-100,500),
breaks = seq(-100,500,100),
guide = "prism_offset"
)
p
image.png
#調(diào)整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
# 調(diào)整x軸和y軸標(biāo)題昼弟，刪除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")
參考鏈接：
