这图怎么画| 多组箱线图+组间/内差异分析

木舟笔记

于 2023-10-07 10:30:06 发布

阅读量457

点赞数

本文链接：https://blog.csdn.net/weixin_45822007/article/details/133663328

版权

box_cover

写在前面

【这图怎么画】系列的图都来自VIP群里同学的提问。推文只是对图片的复现，不代表作者对图片展现形式的认同。欢迎同学们在群里分享有意思的图片。

本期图片

❝
Doi：https://doi.org/10.1111/nph.19115
❞

相似图片：跟着Nat Commun学作图 | Post-hoc图（Extended error bar plot

复现结果

result

❝
可在AI上微调。使用代码时请按需修改，切勿全部套用。
❞

示例数据和代码领取

木舟笔记永久VIP企划

「权益：」

「木舟笔记所有推文示例数据及代码(「在VIP群里」实时更新」）。

data+code
木舟笔记「科研交流群」。

「收费：」

「169￥/人」。可添加微信：mzbj0002 转账（或扫描下方二维码），或直接在文末打赏。木舟笔记「2022VIP」可直接支付「70￥」升级。

❝
点赞、在看 本文，分享至朋友圈集赞30个并保留30分钟，可优惠20￥。
❞

绘图

rm(list = ls())
library(ggplot2) 
library(ggpubr)
library(reshape2) 
library(dplyr) 
# data example
mat <- matrix(runif(40*12, min=0, max=10), nrow=40, ncol=12)
rownames(mat) <- paste0("sample", 1:40)
colnames(mat) <- paste0("gene", 1:12)
mat
gene_group = data.frame(gene = paste0("gene", 1:12),
                        group = rep(paste0("group", 1:6), each = 2))
sample_group = data.frame(sample = paste0("sample", 1:40),
                          group = rep(c('c','w'),20))
# data preprocess
library(reshape2)
mat_long <- melt(mat, varnames = c("sample", "gene"), value.name = "expression")
df1 <- merge(mat_long, gene_group, by = "gene")
df2 <- merge(mat_long, sample_group, by = "sample")

## calculate mean value
df1_mean = df1 %>% 
  group_by(gene) %>% 
  summarise(mean_value=mean(expression))
df1_mean$x = 1:12

df2_mean <- df2 %>%
  group_by(gene, group) %>%
  summarise(mean = mean(expression), sd = sd(expression))
df2_mean$gene <- factor(df2_mean$gene)

# plot
## left boxplot
df1$facet = 'facet1'# for facet title

p1 <- ggplot(df1,aes(gene,expression))+
  geom_boxplot(aes(fill=group,color=group),outlier.shape = 18,size=0.6)+
  stat_boxplot(aes(color=group),geom = "errorbar", width=0.3,size=0.6)+
  geom_segment(df1_mean,
               mapping=aes(x=x-0.5,xend=x+0.5,y=mean_value,yend = mean_value),
               color="white",size=0.5)+
  geom_hline(yintercept = 5, linetype = 2, color = "grey60",linewidth=0.8)+
  ## compare with gene6
  stat_compare_means(label = "p.signif", method = "t.test",
                     ref.group = "gene1",size = 5)+
  coord_flip()+
  scale_fill_manual(values = c("#e3ac6d","#9d7bb8","#6caf83","#d9586e","#3c74bb","#f85b2b"))+
  scale_color_manual(values = c("#e3ac6d","#9d7bb8","#6caf83","#d9586e","#3c74bb","#f85b2b"))+
  theme_bw()+
  theme(legend.position = "none",
        panel.grid = element_blank(),
        axis.text = element_text(color = "black",size=10))+
  labs(y="exp",x='')+
  annotate("rect", xmin = 0, xmax = 4.5, ymin = -Inf, ymax = Inf, alpha = 0.3,fill="#EAA9CF") +
  annotate("rect", xmin = 4.5, xmax = 8.5, ymin = -Inf, ymax = Inf, alpha = 0.3,fill="#D6BDFF") +
  annotate("rect", xmin = 8.5, xmax = 13, ymin = -Inf, ymax = Inf, alpha = 0.3,fill="#99DAFF")+
  facet_grid(~ facet)
p1

## right plot
df2_mean$facet = 'facet2' # for facet title
p2 <- ggplot(df2_mean,aes(x =gene, y = mean,
                          group = group,color = group)) +
  geom_point(position = position_dodge(0.8),shape=18,size=3)+
  geom_errorbar(aes(x =gene, y = mean,ymin = mean- sd, ymax = mean + sd,
                                       group = group,color = group), 
                width = 0,position = position_dodge(0.8),linewidth=0.5)+
  scale_color_manual(values = c("#7fc190","#efb684"))+
  stat_compare_means(df2, mapping = aes(gene,expression,group = group),
                     label = "p.signif", method = "t.test",size = 5)+
  coord_flip()+
  theme_bw()+
  theme(legend.position = c(0.9,0.95),
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.grid = element_blank(),
        axis.text.x = element_text(color = "black",size=10),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())+
  labs(y='exp2',x='',color='')+
  annotate("rect", xmin = 0, xmax = 4.5, ymin = -Inf, ymax = Inf, alpha = 0.3,fill="#EAA9CF") +
  annotate("rect", xmin = 4.5, xmax = 8.5, ymin = -Inf, ymax = Inf, alpha = 0.3,fill="#D6BDFF") +
  annotate("rect", xmin = 8.5, xmax = 13, ymin = -Inf, ymax = Inf, alpha = 0.3,fill="#99DAFF")+
  facet_grid(~ facet)
p2
#merge
library(patchwork)
library(cowplot)
p1+p2+plot_layout(nrow= 1,
                  width = c(3, 2))
ggsave('boxplot.pdf',width = 10,height = 8)