[GCB图表复现] 多组间方差分析热图

2023-12-04 19:34:48 浏览数 (2)

论文

加载R包

代码语言:javascript复制
library(tidyverse)
library(magrittr)
library(multcompView)

数据清洗

代码语言:javascript复制
dff <- read_tsv("phylum.xls") %>%
  pivot_longer(-phylum) %>% left_join(.,read_tsv("group.xls"),by=c("name"="sample")) %>% 
  select(-name)

df <- dff %>% group_by(phylum,group) %>% 
  summarise_all(~sum(.,na.rm=T)) %>% as_tibble()

方差分析

代码语言:javascript复制
p <- split(dff,list(dff$phylum))
aov_data <- data.frame()

str(p)

for (i in 1:22){
  anova <- aov(value ~ group,data=p[i] %>% as.data.frame() %>% 
        set_colnames(c("phylum","value","group")))
  
  Tukey <- TukeyHSD(anova)
  cld <- multcompLetters4(anova,Tukey)
  
  dt <- p[i] %>% as.data.frame() %>% 
    set_colnames(c("phylum","value","group")) %>% 
    group_by(group,phylum) %>%
    summarise(value_mean=mean(value)) %>%
    ungroup() %>% 
    as.data.frame() %>% drop_na()
  
  cld <- as.data.frame.list(cld$`group`)
  dt$Tukey <- cld$Letters
  aov_data  <- rbind(aov_data,dt)

}

数据可视化

代码语言:javascript复制
df %>% left_join(.,aov_data %>% select(1,2,4),by=c("phylum","group")) %>% 
  ggplot(aes(group,phylum,fill=value)) 
  geom_tile() 
  geom_text(aes(label=Tukey)) 
  scale_color_gradientn(colours = rev(RColorBrewer::brewer.pal(3, "RdBu")))  
  scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(3, "RdBu"))) 
  scale_x_discrete(expand = c(0,0),position = "top") 
  scale_y_discrete(expand = c(0,0)) 
  labs(x=NULL,y=NULL) 
  theme(axis.text.x = element_text(angle = 0, hjust = 0, vjust = 0.5,
                                   color = "black", size = 8),
        axis.text.y = element_text(color = "black", size = 8,angle=0),
        axis.ticks = element_blank(),
        panel.spacing.y = unit(0, "cm"),
        plot.background = element_blank(),
        panel.background = element_blank(),
        legend.title = element_blank()) 
  guides(fill = guide_colorbar(direction = "vertical", reverse = F, barwidth = unit(.5, "cm"),
                               barheight = unit(15, "cm")))

0 人点赞