跟着nature communications学绘图(8) ggprism优雅的添加p值

R语言数据分析指南

发布于 2022-09-21 15:10:55

4660

发布于 2022-09-21 15:10:55

❝本节来介绍如何灵活使用「rstatix」来进行数据统计分析，以前都是使用R内置数据来进行展示这次使用论文中的原始数据来进行展示；

加载R包

library(tidyverse)
library(rstatix)
library(ggpubr)
library(ggprism)

数据清洗

df <- read_tsv('F1-a.txt') %>% pivot_longer(-`MUFA-PI / total PI [%]`) %>% 
  group_by(name) %>%
  summarise(value_mean=mean(value),sd=sd(value),se=sd(value)/sqrt(n())) %>% 
  mutate(group=case_when(name=="w/o" ~ "black",name=="TNFα" ~ "grey",
                         name=="I3M" ~ "grey",TRUE ~ "green"))

# A tibble: 10 × 5
   name  value_mean    sd    se group
   <fct>      <dbl> <dbl> <dbl> <chr>
 1 CHX         19.6  2.45 1.22  green
 2 ETO         32.7  2.46 1.23  green
 3 I3M         50.1  3.51 1.76  grey 
 4 MC          39.8  4.85 2.43  green
 5 Serum       22.4  3.36 1.68  green
 6 STS         21.7  1.16 0.578 green
 7 TNFα        52.7  5.97 2.98  grey 
 8 TPG         35.1 11.0  5.49  green
 9 VAL         31.1  5.60 2.80  green
10 w/o         45.1  3.29 1.65  black

定义因子

df$name <- factor(df$name,levels=read_tsv('F1-a.txt') %>% 
                    pivot_longer(-`MUFA-PI / total PI [%]`) %>%
                    select(name) %>% distinct() %>% pull())

统计分析

❝此处还是通过整合数据得到最高点位置信息，如果使用代码生成则会存在问题 ❞

stat.test <- read_tsv('F1-a.txt') %>% pivot_longer(-`MUFA-PI / total PI [%]`) %>% 
  t_test(data =., value ~ name, ref.group  =  "w/o") %>% 
  mutate(p.adj.signif = replace_na(p.adj.signif,""),across("p.adj.signif",str_replace,"ns","")) %>% 
  select(group1,group2,p.adj,p.adj.signif) %>% 
  left_join(.,df,by=c("group2"="name")) %>% 
  mutate(y.position=value_mean+se+0.3)

定义主题

theme_niwot <- function(){
  theme_test() +
    theme(axis.title = element_blank(),
          axis.ticks.x=element_blank(),
          panel.grid.major.y = element_line(color = "#DAE1E7"),
          panel.grid.major.x = element_blank(),
          plot.margin = unit(rep(0.2,4),"cm"),
          axis.text = element_text(size = 10, color = "#22292F"),
          axis.text.y = element_text(margin = margin(r = 5)),
          axis.text.x = element_text(margin = margin(t = 5)),
          legend.position = "non")
}

数据可视化

df %>% ggplot(.,aes(name,value_mean))+
  geom_errorbar(aes(ymax = value_mean + se, ymin = value_mean - se),width = 0.2,color = "grey30")+
  geom_col(width=0.5,aes(fill=group),color="grey50")+
  add_pvalue(stat.test,label = "p.adj.signif",label.size=6,
             coord.flip = TRUE, remove.bracket = TRUE)+
  scale_y_continuous(expand=c(0,0),limits = c(0,57)) +
  theme_niwot()+
  scale_fill_brewer(palette="Blues")