R语言计算多样性指数的平均值和方差

在16s RNA高通量测序后，我们经常要对样本进行Alpha多样性指数分析，所以经常要计算Alpha多样性指数的平均值和方差，如果在你的分析中需要用到下图中这种表示方法，你可以参考该文章，否则跳过。

示例表格

1.Alpha 多样性指数原始数据的准备

alpha文本文件 是一组经过测序后经过计算得到的alpha多样性指数，包括ACE、Chao1、Simpson、Shannon等多种指数，这里就不做详细讲解。

group文本文件是一个样本的分组文件，包括样本名、分组1、分组2.

图1-1 alpha文本文件

图1-2 group文本文件

2 R语言计算多样性指数的平均值和方差

2.1 数据预处理

下面我们将使用R语言来计算多样性指数的平均值和方差，当然你也可以使用excel来做，也很方便。

#设置工作目录####
setwd("D:/R_wenji/06-微信公众号/21_07_03R语言计算多样性指数的平均值和方差")
#导入Alpha.txt
data<-read.table("Alpha.txt",header=T,sep="\t",row.names=1)
#  提取ACE、Chao1、Simpson、Shannon四种指数
data <- data[,1:4] 
#导入group分组文件
group  <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)

图2.1-1 导入Alpha.txt

图2.1-2 提取ACE、Chao1、Simpson、Shannon四种指数

图2.1-3 导入group分组文件

#添加样本名
data$sample<- factor(rownames(data), levels = rev(rownames(data)))
#合并两个表格
data <- merge(data, group, by = 'sample')
#将site转化为因子变量
data$site <- factor(data$site)
str(data) #查看数据类型

图2.1-4 添加样本名

图2.1-5 合并两个表格

图2.1-6 查看数据类型

2.2 R语言计算多样性指数的平均值和方差

#计算ACE平均值
ACE_mean <- aggregate(data$ACE, by = list(data$site), FUN = mean)
names(ACE_mean) <- c('sample','ACE_mean' ) #修改列名
#保留两位小数
ACE_mean$ACE_mean <- sprintf("%0.2f", ACE_mean$ACE_mean)

#计算ACE标准差
ACE_sd <- aggregate(data$ACE, by = list(data$site), FUN = sd)
names( ACE_sd)<-  c('sample','ACE_sd' ) #修改列名
ACE_sd$ACE_sd <- sprintf("%0.2f", ACE_sd$ACE_sd)
#合并两个表格
ACE <- merge(ACE_mean, ACE_sd, by = 'sample')

图2.2-1 ACE

#加载tidyr包，没安装需要先安装该包
library(tidyr) 
#添加列mean±sd
#ACE1 <- tidyr::unite(ACE,"ACE_mean ± ACE_sd", ACE_mean, ACE_sd,sep = "±")
ACE <- unite(ACE, "ACE_mean ± ACE_sd", ACE_mean, ACE_sd, sep = "±", remove = FALSE)

图2.2-2 添加ACE_mean ± ACE_sd列

2.3 其他多样性指数方法和ACE计算一样

#计算chao1标准差和平均值
chao1_mean <- aggregate(data$Chao1, by = list(data$site), FUN = mean)
names(chao1_mean) <- c('sample','chao1_mean' ) #修改列名
chao1_mean$chao1_mean <- sprintf("%0.2f", chao1_mean$chao1_mean)
chao1_sd <- aggregate(data$Chao1, by = list(data$site), FUN = sd)
names( chao1_sd)<-  c('sample','chao1_sd' ) #修改列名
chao1_sd$chao1_sd <- sprintf("%0.2f", chao1_sd$chao1_sd)
#合并两个表格
chao1 <- merge(chao1_mean, chao1_sd, by = 'sample')
chao1 <- unite(chao1, "chao1_mean ± chao1_sd", chao1_mean, chao1_sd, sep = "±", remove = FALSE)

#计算Simpson标准差和平均值
Simpson_mean <- aggregate(data$Simpson, by = list(data$site), FUN = mean)
names(Simpson_mean) <- c('sample','Simpson_mean' ) #修改列名
Simpson_mean$Simpson_mean <- sprintf("%0.2f", Simpson_mean$Simpson_mean)

Simpson_sd <- aggregate(data$Simpson, by = list(data$site), FUN = sd)
names(Simpson_sd)<-  c('sample','Simpson_sd' ) #修改列名
Simpson_sd$Simpson_sd <- sprintf("%0.2f", Simpson_sd$Simpson_sd)
#合并两个表格
Simpson <- merge(Simpson_mean,Simpson_sd, by = 'sample')

Simpson <- unite(Simpson, "Simpson_mean ± Simpson_sd", Simpson_mean, Simpson_sd, sep = "±", remove = FALSE)

#计算Shannon标准差和平均值
Shannon_mean <- aggregate(data$Shannon, by = list(data$site), FUN = mean)
names(Shannon_mean) <- c('sample','Shannon_mean' ) #修改列名
Shannon_mean$Shannon_mean <- sprintf("%0.2f", Shannon_mean$Shannon_mean)

Shannon_sd <- aggregate(data$Shannon, by = list(data$site), FUN = sd)
names(Shannon_sd)<-  c('sample','Shannon_sd' ) #修改列名
Shannon_sd$Shannon_sd <- sprintf("%0.2f", Shannon_sd$Shannon_sd)

#合并两个表格
Shannon<- merge(Shannon_mean,Shannon_sd, by = 'sample')
Shannon <- unite(Shannon, "Shannon_mean ± Shannon_sd", Shannon_mean, Shannon_sd, sep = "±", remove = FALSE)

计算完四个指数后，我们需要将他们合并在一张表格

#合并所有指数的结果
Alpha1<- merge(ACE,chao1, by = 'sample')
Alpha1<- merge(Alpha1,Simpson, by = 'sample')
Alpha1<- merge(Alpha1,Shannon, by = 'sample')

#查看其数据类型
str(Alpha1)

图2.3

#将数据导出
write.table (Alpha1, file ="Alpha_1.csv",sep =",", quote =FALSE) #将数据导出

我们可以使用代码将结果导出，然后在word中修改为三线表，好了今天的内容就是这些，我感觉比excel方便一些，你觉得呢？

感谢你的阅读，谢谢你的支持，未来可期，加油。