箱线图进行方差分析并添加显著性标记

❝本节来介绍如何「在计算多样性指数的基础上来进行显著性标记」；

加载R包

代码语言：javascript复制

library(tidyverse)
library(vegan)
library(magrittr)
library(multcompView)

导入数据

代码语言：javascript复制

alpha <- read.delim("otu_taxa_table-2.xls",sep="t",row.names = 1) %>% 
  t() %>% as.data.frame()

group <- read_tsv("group.xls") %>% set_colnames(c("sample","group"))

定义函数计算多样性指数

代码语言：javascript复制

alpha_diversity <- function(x,y) {
  Shannon <- diversity(x, index = 'shannon')
  Simpson <- diversity(x, index = 'simpson')  
  observed_species <- specnumber(x)
  Chao1 <- estimateR(x)[2,]
  ACE <- estimateR(x)[4,]
  pielou <- diversity(x,index = "shannon")/log(specnumber(x),exp(1))
  
  result <- data.frame(Shannon,Simpson,observed_species,Chao1,ACE,pielou) %>% 
    rownames_to_column("sample") %>% 
    left_join(.,y,by="sample")
    
  return(result)
  
}

数据整理

代码语言：javascript复制

df <- alpha_diversity(alpha,group) %>% select(-sample,-observed_species,-Simpson) %>% 
  pivot_longer(-group)

定义颜色

代码语言：javascript复制

col <- c("#1F78B4","#33A02C","#FB9A99","#E31A1C","#FDBF6F","#B2DF8A",
         "#A6CEE3","#BA7A70","#9D4E3F","#829BAB")

❝上面这些基本是上一篇文档的内容为了文档结构的完整，将其放置于此；那么接下来就是本文的重点内容多组之间进行方差分析添加显著性标记 ❞

方差分析

代码语言：javascript复制

p <- split(df,list(df$name))
aov_data <- data.frame()

str(p)

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

构建显著性标记数据集

代码语言：javascript复制

df2 <- df %>% arrange(name) %>% left_join(.,aov_data,by=c("group","name"))

text <- df2 %>% group_by(group,name) %>% summarise(max(value)) %>% arrange(name) %>% ungroup() %>% 
  set_colnames(c("group","name","value")) %>% 
  left_join(.,df2 %>% select(1,2,6),by=c("group","name")) %>% distinct() %>% 
  mutate(value=case_when(name =="ACE" ~ value 90,
                         name =="Chao1" ~ value 90,
                         name =="pielou" ~ value  0.008,
                         name =="Shannon" ~ value 0.065))

❝由于循环构建的为条形图的数据，但显著性标记是不区分图形的因此在此通过上面的代码构建箱线图的数据，由于还存在离群值因此做了过多的处理，各位观众老爷细细品味 ❞

定义绘图函数

代码语言：javascript复制

make_plot <- function(data,x,y,z){
  ggplot(data,aes(x={{x}},y={{y}},fill={{x}})) 
    stat_boxplot(geom="errorbar",position=position_dodge(width=0.2),width=0.2) 
    geom_boxplot(position=position_dodge(width =0.2),width=0.5,outlier.shape = NA) 
    scale_fill_manual(values={{z}}) 
    facet_wrap(.~name,scales = "free") 
    theme_bw() 
    theme(panel.spacing.x = unit(0.2,"cm"),
          panel.spacing.y = unit(0.1, "cm"),
          axis.title = element_blank(),
          strip.text.x = element_text(size=12,color="black"),
          axis.text = element_text(color="black"),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          legend.position = "non",
          plot.margin=unit(c(0.3,0.3,0.3,0.3),units=,"cm"))
}

数据可视化

代码语言：javascript复制

make_plot(df,group,value,col) 
  geom_text(data=text,aes(label=Tukey,y=value))

数据可视化

0 人点赞