R语言实现Alpha多样性指数的计算
上次我们已经使用R语言来对OTU表的抽平分析，那么我们如何使用抽平后的OTU表来重新计算Alpha多样性呢？接下来我们就来学习一下。
不过你可能会说，这个不是测序公司都计算好了吗，而且还可以用qiime软件，为啥还要使用R，因为有些测序公司并没有帮你抽平，再一个我使用R语言比较多，如果你会使用其他软件当然更好。欢迎你与我分享一下。
1 所需的数据类型
这里我们需要使用到两个数据集：一个是经过抽平分析处理后的otu表（如果不会抽平分析的可以查看该文章
；
另外一个数据集是使用各OTU代表序列构建的进化树文件“otu_tree.tre”。
抽平分析后得到的otu表：（当然你也可以选择不抽平）
图1-1 otu表
图1-1 otu_tree.tre文件
计算每一种Alpha多样性指数都会用到otu表，但是对于otu_tree.tre文件，只用于计算谱系多样性。
2 使用R语言计算常用的Alpha多样性指数
我们接下来会使用到两个包，一个是vegan包，另一个是picante包。如果没安装这两个包，需要提前安装好。
vegan包可以用来计算多种Alpha多样性指数，例如这次我们要学习计算的物种丰富度（Richness）、Chao 1指数、ACE指数、Shannon指数、Simpson指数等。
谱系多样性（即PD_whole_tree）需要使用picante包，该多样性除了物种丰富度数据外还需要进化树文件。
2.1 加载包以及数据集
#设置工作目录
setwd("D:/R_wenji/06-微信公众号/21_07_05")
#需要加载vegan包和picante包，没有安装需要先安装
library(vegan)
library(picante)
#读入抽平后的otu表
otu <- read.delim('otu.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#将otu数据转置
otu <- t(otu)
#加载进化树文件
tree <- read.tree('otu_tree.tre')
图2.1-1 读入抽平后的otu表
2.2 计算-物种丰富度 Richness 指数，又称observed species 指数
#计算方法一：
observed_species <- rowSums(otu > 0)
#计算方法二：
observed_species <- estimateR(otu)[1, ]
#输出查看结果
observed_species
图2.2 物种丰富度 Richness 指数
2.3 计算 Chao 1指数
Chao1 <- estimateR(otu)[2, ]
Chao1
图2.3 Chao 1指数
2.4 计算ACE 指数
ACE <- estimateR(otu)[4, ]
ACE
图2.4 ACE指数
2.5 计算Shannon指数
#Shannon 指数,通常使用2、e作为底数
#以e作为底数表示方法
Shannon <- diversity(otu, index = 'shannon', base = exp(1))
#以2作为底数表示方法
Shannon <- diversity(otu, index = 'shannon', base = 2)
#输出Shannon_index结果
Shannon
图2.5 Shannon指数
2.6 计算Simpson指数
#Simpson指数分为经典 Simpson 指数和Gini-Simpson 指数，不过平时常用的 Simpson 指数即为 Gini-Simpson 指数
#Gini-Simpson 指数代码
Gini_simpson <- diversity(otu, index = 'simpson')
#经常使用
Gini_simpson
#经典 Simpson 指数
simpson_index <- 1 - Gini_simpson
图2.6 Simpson指数
2.7 计算goods_coverage 指数
goods_coverage <- 1 - rowSums(otu == 1) / rowSums(otu)
goods_coverage
图2.7 goods_coverage 指数
2.8 计算谱系多样性（PD）
#除了otu文件，需要指定进化树文件
PD_whole_tree <- pd(otu, tree, include.root = FALSE)[1]
PD_whole_tree
图2.8 谱系多样性（PD）
单个计算Alpha多样性指数的方法已经讲完了，那么我该如何使用会比较方便呢？
那么请使用下面这个自定义函数，函数不够完美，你有需求可以自己修改。
3 使用自定义alpha_diversity函数来快速计算多种Alpha多样性指数
首先我们需要定义alpha_diversity函数：
library(vegan)
library(picante)
alpha_diversity <- function(x, tree = NULL) {
observed_species <- estimateR(x)[1, ]
Chao1 <- estimateR(x)[2, ]
ACE <- estimateR(x)[4, ]
Shannon <- diversity(x, index = 'shannon',base = 2)
Simpson <- diversity(x, index = 'simpson') #注意，这里是Gini-Simpson 指数
goods_Coverage <- 1 - rowSums(x == 1) / rowSums(x)
#保留四位小数
Shannon <- sprintf("%0.4f", Shannon)
Simpson <- sprintf("%0.4f", Simpson)
goods_Coverage <- sprintf("%0.4f", goods_Coverage)
result <- data.frame(observed_species, ACE,Chao1, Shannon, Simpson, goods_Coverage)
if (!is.null(tree)) {
PD_whole_tree <- pd(x, tree, include.root = FALSE)[1]
names(PD_whole_tree) <- 'PD_whole_tree'
result <- cbind(result, PD_whole_tree)
result <- data.frame(observed_species, ACE,Chao1, Shannon, Simpson,
PD_whole_tree ,goods_Coverage)
}
result
}
alpha_diversity函数定义好了，我们就可以导入数据进行计算了
#加载OTU 表
otu <- read.delim('otu.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- t(otu)
#加载进化树文件
tree <- read.tree('otu_tree.tre')
#如果不需要计算谱系多样性
alpha <- alpha_diversity (otu)
#需要计算谱系多样性时，需要指定进化树文件
alpha1 <- alpha_diversity (otu, tree)
图3-1 不包括谱系多样性的Alpha多样性指数
图3-2 包括谱系多样性的Alpha多样性指数
#将结果输出，保存在本地
write.csv(alpha, 'alpha_diversity.csv', quote = FALSE)
write.csv(alpha1, 'alpha_diversity1.csv', quote = FALSE)
图3-3 你可以在文件夹中看到它们
是不是非常简单，如果你需要获取示例数据及代码，可以给我留言，如果你觉得对你有帮助，记得点个赞。
让我们一起加油吧。
