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 你可以在文件夹中看到它们

是不是非常简单，如果你需要获取示例数据及代码，可以给我留言，如果你觉得对你有帮助，记得点个赞。

让我们一起加油吧。