vlambda博客
学习文章列表

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

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_coveragePD_whole_tree8种常见的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)

alpha_all”如下所示,对于每个样本,计算了多种Alpha多样性指数供后续分析使用。

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

 


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