提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  其他

当前栏目

跟着Nature Communications学作图：R语言ggplot2做柱形图并添加误差线和显著性P值

2023-03-14 22:48:00 时间

论文是

A giant NLR gene confers broad-spectrum resistance to Phytophthora sojae in soybean

image.png

论文里公布了大部分柱形图和箱线图的原始数据，今天的推文试着用论文中的数据模仿一下论文中的 Figure 2b c

image.png

Figure 2b 的数据

image.png

类似的图之前录制过视频进行介绍，如果习惯看视频的话可以关注下我的B站账号 小明的数据分析笔记本

image.png

首先读取数据 

dfb<-read.csv("figure2b.csv",header=F)
dfb

宽格式数据转换为长格式 

dfb %>% 
  pivot_longer(!V1) %>% 
  select(V1,value) %>% 
  na.omit() -> dfb.1

最基本的柱形图 

library(ggplot2)

ggplot(data=dfb.1,aes(x=V1,y=value))+
  stat_summary(geom = "bar",
               fun = mean,
               fill="#c6c3c3")

image.png

添加误差线

这里误差线采用的是mean+-sem

library(ggplot2)

ebtop<-function(x){
  return(mean(x)+sd(x)/sqrt(length(x)))
}
ebbottom<-function(x){
  return(mean(x)-sd(x)/sqrt(length(x)))
}

ggplot(data=dfb.1,aes(x=V1,y=value))+
  stat_summary(geom = "bar",
               fun = mean,
               fill="#c6c3c3")+
  stat_summary(geom = "errorbar",
               fun.min = ebbottom,
               fun.max = ebtop,
               width=0.2)

image.png

添加图上的散点 

library(ggplot2)

ebtop<-function(x){
  return(mean(x)+sd(x))
}
ebbottom<-function(x){
  return(mean(x)-sd(x))
}

ggplot(data=dfb.1,aes(x=V1,y=value))+
  stat_summary(geom = "bar",
               fun = mean,
               fill="#c6c3c3")+
  stat_summary(geom = "errorbar",
               fun.min = ebbottom,
               fun.max = ebtop,
               width=0.2)+
  geom_jitter(width = 0.3)

image.png

添加显著性p值 

ggplot(data=dfb.1,aes(x=V1,y=value))+
  stat_summary(geom = "bar",
               fun = mean,
               fill="#c6c3c3")+
  stat_summary(geom = "errorbar",
               fun.min = ebbottom,
               fun.max = ebtop,
               width=0.2)+
  geom_jitter(width = 0.3)+
  geom_signif(comparisons = list(c("Control","F5"),
                                 c("Control","pAtUbi3:CDS-Rps11-1"),
                                 c("Control","pAtUbi3:CDS-Rps11-2")),
              test = t.test,
              test.args = list(var.equal=T,
                               alternative="two.side"),
              y_position = c(1.1,1.3,1.5),
              #annotations = c(""),
              parse = T)

image.png

如何在geom_signif()函数里调整P值的文字格式暂时想不到办法了，使用annotate()函数吧

ggplot(data=dfb.1,aes(x=V1,y=value))+
  stat_summary(geom = "bar",
               fun = mean,
               fill="#c6c3c3")+
  stat_summary(geom = "errorbar",
               fun.min = ebbottom,
               fun.max = ebtop,
               width=0.2)+
  geom_jitter(width = 0.3)+
  geom_signif(comparisons = list(c("Control","F5"),
                                 c("Control","pAtUbi3:CDS-Rps11-1"),
                                 c("Control","pAtUbi3:CDS-Rps11-2")),
              test = t.test,
              test.args = list(var.equal=T,
                               alternative="two.side"),
              y_position = c(1.1,1.3,1.5),
              annotations = c(""),
              parse = T)+
  annotate(geom = "text",
           x=1.5,y=1.15,
           label=expression(italic(P)~'='~1.83%*%10^-6))+
  annotate(geom = "text",
           x=2,y=1.35,
           label=expression(italic(P)~'='~2.71%*%10^-5))+
  annotate(geom = "text",
           x=2.5,y=1.55,
           label=expression(italic(P)~'='~5.75%*%10^-8))

image.png

这里遇到的警告信息暂时搞不懂是什么意思了

image.png

接下来是细节的调整 

ggplot(data=dfb.1,aes(x=V1,y=value))+
  stat_summary(geom = "bar",
               fun = mean,
               fill="#c6c3c3")+
  stat_summary(geom = "errorbar",
               fun.min = ebbottom,
               fun.max = ebtop,
               width=0.2)+
  geom_jitter(width = 0.3)+
  geom_signif(comparisons = list(c("Control","F5"),
                                 c("Control","pAtUbi3:CDS-Rps11-1"),
                                 c("Control","pAtUbi3:CDS-Rps11-2")),
              test = t.test,
              test.args = list(var.equal=T,
                               alternative="two.side"),
              y_position = c(1.1,1.3,1.5),
              annotations = c(""),
              parse = T)+
  annotate(geom = "text",
           x=1.5,y=1.15,
           label=expression(italic(P)~'='~1.83%*%10^-6))+
  annotate(geom = "text",
           x=2,y=1.35,
           label=expression(italic(P)~'='~2.71%*%10^-5))+
  annotate(geom = "text",
           x=2.5,y=1.55,
           label=expression(italic(P)~'='~5.75%*%10^-8))+
  scale_y_continuous(expand = c(0,0),
                     limits = c(0,1.6),
                     breaks = seq(0,1,0.2))+
  theme_minimal()+
  theme(panel.grid = element_blank(),
        axis.line.y = element_line(),
        axis.ticks.y = element_line(),
        axis.title.y = element_text(hjust=0.25,
                                    size=15),
        axis.text.x = element_text(angle = 30,
                                   hjust = 1,
                                   size=10))+
  guides(y=guide_axis_truncated(trunc_lower = 0,
                               trunc_upper = 1))+
  labs(x=NULL,y="Survival Rate")

image.png

Figure 2c 的数据也有，大家可以试着用以上代码试着复现一下figure2c的数据

相关文章