跟着Nature Genetics 学画图：R语言ggplot2画箱线图(boxplot)展示D statistic

用户7010445

发布于 2021-05-07 10:54:13

2.5K0

发布于 2021-05-07 10:54:13

首先推荐一个分享R语言干货的优质公众号 R语言统计与绘图

简介：R语言统计与绘图公众号目前致力于分享医学统计与R绘图知识，手把手教你使用R语言绘制基线特征表、KM生存曲线、森林图、ROC曲线等。每天一篇精彩R语言推文教程，手把手带你入门R语言绘图。
今天推文重复的图来自于论文 Whole-genome resequencing of 445 Lactuca accessions reveals the domestication history of cultivated lettuce

image.png

这篇论文的数据是公开的，代码也公开了一部分，那我们就可以按照他的提供的数据来试着复原一些论文中的图了。本来已经重复到了论文中Fig1的2c，但是试着做局部放大的时候遇到了一些问题，暂时还搞不定，看了其他的图实现起来还有一定的难度。所以先挑一个相对比较简单的箱线图来模仿一下吧。对应的是论文中的Figure 2d

image.png

数据是论文附件的 Source Data Fig.2

image.png

首先是读入数据

df<-readxl::read_excel("NG/41588_2021_831_MOESM5_ESM.xlsx",
                      sheet="Fig2d")

论文中提供的是宽格式数据，如果使用ggplot2作图需要转换成长格式，这里本来想尝试一下tidyr包中的pivot_longer()函数了，帮助文档没有看明白。没有搞定，还是直接使用reshape2中的melt()函数吧

library(dplyr)
df %>% 
  mutate(new_col=paste(Group1,Group2,sep="_")) %>% 
  select(-c("Group1","Group2","Group3","Outgroup")) %>% 
  #reshape2::melt(var.ids=c("Group1")) %>% 
  #arrange(Group1,Group2) %>% 
  reshape2::melt(var.ids="new_col") -> df1

head(df1)

ggplot2 作图

library(ggplot2)
library(stringr)
library(ggprism)
x_level<-paste(df$Group1,df$Group2,sep="_")
x_level
df1$group<-str_sub(df1$new_col,5,7)
df1$new_col<-factor(df1$new_col,
                    levels = x_level)
ggplot(df1,aes(x=new_col,y=value))+
  stat_boxplot(geom = "errorbar",width=0.2)+
  geom_boxplot(outlier.shape = 1,
               aes(fill=group),
               show.legend = F)+
  scale_fill_manual(values = c("#e64b35",
                               "#4daf4a",
                               "#4dbbd5",
                               "#cab2d6",
                               "#b2df8a"))+
  scale_x_discrete(labels=str_sub(x_level,1,3),
                   guide = "prism_offset")+
  scale_y_continuous(limits = c(-0.021,0.085),
                     breaks = seq(-0.02,0.08,by=0.02))+
  theme_prism(axis_text_angle = 90,
              base_line_size = 0.1,
              base_fontface = "plain",
              base_family = "serif")+
  labs(x=NULL,
       y=expression(paste(italic("D")," statistic")))+
  theme(plot.margin = unit(c(0.2,0.2,2,0.2),'cm'))+
  
  geom_segment(x=1,xend=5,y=-0.04,yend=-0.04)+
  annotate("text",x=3,y=-0.02,label="CAU",vjust=10)+
  geom_segment(x=6,xend=9,y=-0.04,yend=-0.04)+
  annotate("text",x=7.5,y=-0.02,label="SEU",vjust=10)+
  geom_segment(x=10,xend=12,y=-0.04,yend=-0.04)+
  annotate("text",x=11.5,y=-0.02,label="WEU",vjust=10)+
  geom_segment(x=13,xend=14,y=-0.04,yend=-0.04)+
  annotate("text",x=13.5,y=-0.02,label="EEU",vjust=10)+
  geom_segment(x=14.5,xend=15.5,y=-0.04,yend=-0.04)+
  annotate("text",x=15,y=-0.02,label="WAS",vjust=10)+
  annotate("text",x=1,y=-0.02,label="P1",
           hjust=2,vjust=5)+
  annotate("text",x=1,y=-0.02,label="P2",
           hjust=2,vjust=10)+
  coord_cartesian(clip = "off")