论文
加载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")))