跟着Nature Communications学作图：R语言gg

论文是

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的数据

如果需要推文的示例数据和代码的话，直接给推文打赏1元。如果没有收到我回复的数据代码下载链接，可以加我的微信催我 ，我的微信是 mingyan24

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！