欢迎关注R语言数据分析指南

本节来介绍如何在计算多样性指数的基础上来进行显著性标记，可在文末找到获取数据的方式

加载R包

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

导入数据

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"))

定义函数计算多样性指数

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)
  
}

数据整理

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

定义颜色

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

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

方差分析

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)
}

构建显著性标记数据集

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))

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

定义绘图函数

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"))
}

数据可视化

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

数据获取

可以看到过程还是比较繁琐的没有直接调用R包来的方便，但是通过此文你一定有所收获，需要获取数据的可以关注我的公众号R语言数据分析指南会找到数据获取方式