R语言计算群落多样性分析中常见Alpha多样性指数

vlambda
2020-05-27

R语言计算群落多样性分析中常见Alpha多样性指数

R语言计算Alpha多样性指数

上篇已经对群落多样性分析中的几种常见的概念作了普及。

本篇以某16S高通量测序所得细菌群落数据为例，简介使用R计算几种在微生物群落分析中常用的Alpha多样性指数。

本篇示例文件、R脚本等，已上传至百度盘，提取码xim2。

https://pan.baidu.com/s/1fUiG1WMK6T8cpAQg4trXDg

网盘文件“otu_table.txt”为某16S扩增子测序所得OTU丰度表格，可以理解为物种丰度表。表中每一列为一个样本，每一行为一种OTU，交叉区域为每种OTU在各样本中的丰度。

R语言计算群落多样性分析中常见Alpha多样性指数

另一文件“otu_tree.tre”为使用各OTU代表序列构建的进化树文件（即OTU水平的16S进化树，已存储为nwk格式）。纯文本模式打开查看该进化树文件，内容长这样。

R语言计算群落多样性分析中常见Alpha多样性指数

接下来使用R，基于上述两个文件，计算微生物群落的多种Alpha多样性指数。每一种Alpha多样性指数都会用到OTU丰度表；对于进化树文件，将用于计算一种特殊的多样性指数，谱系多样性。

R语言计算常见的Alpha多样性指数

加载R包，读入OTU丰度表，计算常见Alpha多样性指数，概念。

vegan包是生态群落分析的常用包，该包提供了计算多种Alpha多样性指数的命令，加载vegan包使用。常见的仅基于物种丰度数据即可得到的多样性指数，例如物种丰富度（Richness）、Shannon指数、Simpson指数、均匀度、Chao1指数、ACE指数等，均可基于vegan包中的命令计算得到。

谱系多样性（即PD_whole_tree）无法基于vegan包计算出，因为该指数比较特殊，除了物种丰富度外还需给定进化树文件。使用picante包可计算谱系多样性。

library(vegan)
library(picante)
#注：picante 包加载时默认同时加载 vegan，如果加载了它，可省略“library(vegan)”这一步
 
#读入物种数据
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- t(otu)
 
##物种丰富度 Richness 指数
richness <- rowSums(otu > 0)
#或
richness <- estimateR(otu)[1, ]
 
##Shannon（以下 Shannon 公式的对数底数均设使用 e，在 R 中即表示为 exp(1)）
#Shannon 指数
shannon_index <- diversity(otu, index = 'shannon', base = exp(1))
 
#Shannon 多样性
shannon_diversity <- exp(1)^shannon_index
 
#Shannon 均匀度（Pielou 均匀度）
pielou <- shannon_index / log(richness, exp(1))
 
##Simpson
#Gini-Simpson 指数（我们平时常用的 Simpson 指数即为 Gini-Simpson 指数）
gini_simpson_index <- diversity(otu, index = 'simpson')
 
#经典 Simpson 指数（使用频率比较低）
simpson_index <- 1 - gini_simpson_index
 
#Invsimpson 指数（Gini-Simpson 的倒数）
invsimpson_index <- 1 / gini_simpson_index
#或
invsimpson_index <- diversity(otu, index = 'invsimpson')
 
#Simpson 多样性
simpson_diversity <- 1 / (1 - gini_simpson_index)
 
#Simpson 均匀度（equitability 均匀度）
equitability <- 1 / (richness * (1 - gini_simpson_index))
 
##Chao1 & ACE
#Chao1 指数
chao1 <- estimateR(otu)[2, ]
 
#ACE 指数
ace <- estimateR(otu)[4, ]
 
##goods_coverage 指数
goods_coverage <- 1 - rowSums(otu == 1) / rowSums(otu)
 
##谱系多样性（与上述相比，还需指定进化树文件）
#测试时发现有根树和无根树的PD_whole_tree计算结果是一样的，但是无根树的计算会更快
tree <- read.tree('otu_tree.tre')
pd_whole_tree <- pd(otu, tree, include.root = FALSE)

熟悉了计算各Alpha多样性指数的命令后，接下来我们考虑将它们整合在一起，方便一步得到多个Alpha多样性指数结果。

整合上述命令，定义一个统一的函数alpha()，包含丰富度、Shannon熵指数、Simpson指数（Gini-Simpson指数）、Pielou均匀度、Chao1指数、ACE指数、goods_coverage、PD_whole_tree共8种常见的Alpha多样性指数。

#定义函数
library(picante)       #picante 包加载时默认同时加载 vegan
 
alpha <- function(x, tree = NULL, base = exp(1)) {
        est <- estimateR(x)
        Richness <- est[1, ]
        Chao1 <- est[2, ]
        ACE <- est[4, ]
        Shannon <- diversity(x, index = 'shannon', base = base)
        Simpson <- diversity(x, index = 'simpson')    #Gini-Simpson 指数
        Pielou <- Shannon / log(Richness, base)
        goods_coverage <- 1 - rowSums(x == 1) / rowSums(x)
        
        result <- data.frame(Richness, Shannon, Simpson, Pielou, Chao1, ACE, 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
}
 
#现在直接使用定义好的命令 alpha()，一步得到多种 Alpha 多样性指数
#加载 OTU 丰度表和进化树文件
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- t(otu)
tree <- read.tree('otu_tree.tre')
 
#不包含谱系多样性，无需指定进化树；Shannon 公式的 log 底数我们使用 2
alpha_all <- alpha(otu, base = 2)
#包含谱系多样性时，指定进化树文件；Shannon 公式的 log 底数我们使用 2
alpha_all <- alpha(otu, tree, base = 2)
 
#输出保存在本地
write.csv(alpha_all, 'alpha.csv', quote = FALSE)