alpha多样性指标计算

#以下内容引自https://mp.weixin.qq.com/s?__biz=MzIxNzc1Mzk3NQ==&mid=2247484134&idx=1&sn=345bcbd599b01059db167ee3c423930a&chksm=97f5b2fea0823be83a558ef30a6d54e870a030a031e3fe3bfa540320015c2c6717b6b6a1f3b1&cur_album_id=1886058347326423044&scene=189#wechat_redirect

Whittaker（1972）提出了以下三种多样性。
Alpha多样性（α多样性），指某个群落或生境内部的种的多样性。
Beta多样性（β多样性），即在一个梯度上从一个生境到另一个生境所发生的种的多样性变化的速率和范围。它是研究群落之间的种多度关系。
Gamma多样性（γ多样性），即在一个地理区域内（例如一个岛屿）一系列生境中种的多样性。它通过这些生境的α多样性和生境之间的β多样性的研究范围结合起来表示。

（图片来源，Jurasinski and Beierkuhnlein，2009）
对于Alpha多样性而言，常通过一系列的指数来表征。通过比较这些Alpha多样性指数的高低，判断群落物种多样性水平。
一、物种丰富度
物种丰富度指数（Species richness）为群落中丰度大于0的物种数之和，值越大表明群落中物种种类越丰富。下式中，S，物种丰富度指数；n，个体数（丰度）大于0的物种类型总数。
物种丰富度指数在计算中对所有存在的物种（无论优势物种或稀有物种）等权重看待，只关注物种存在与否，与它们的相对丰度无关。此外，丰富度指数对抽样深度所造成的差异也是非常敏感的。
在微生物16S/18S/ITS等测序分析中，常见observed species指数，即为物种丰富度指数。
二、香农指数（Shannon）
香农指数（Shannon index）或称香农熵指数（Shannon entropy index）、香农-威纳指数（Shannon-Wiener index），同时考虑了物种丰富度以及均匀度（Shannon，1948a，b）。它由信息论（information theory）延伸而来，反映了我们能够预测在群落中随机选择的个体属于哪些物种的不确定性。如果群落仅由单一物种组成（种群），那么我们确信随机选择的个体必定为那个唯一的物种，此时不确定性就为零；否则，我们将无法得知随机被选择的个体究竟属于什么物种，并且不确定性也会随着群落物种种类数的增多而增加。但是，如果群落中存在一种或少数几种物种占据了优势地位（与其它物种相比，它们在丰度上具有明显的优势），那么不确定性就不会那么高，因为我们随机选择的个体很有可能就是这些优势物种。
三、辛普森指数（Simpson）
辛普森指数（Simpson index）同样考虑了物种丰富度以及均匀度，但与Shannon指数相比，它更受均匀度的影响（Simpson，1949)）。经典Simpson指数代表了在群落中两个随机选择的个体属于同一物种的概率，当群落物种丰富度增加时，这种概率降低，即Simpson指数随着物种丰富度的增加而降低。由于经典Simpson指数与物种丰富度相反的趋势不直观，如今常用演变而来的Gini-Simpson指数表示Simpson指数，即用1减去经典Simpson指数数值后得到，此时Simpson指数随着丰富度的增加而增加（二者保持一致的趋势）。
通常情况下为保持物种丰度指数与Simpson指数的趋势一致，直接将“Gini-Simpson”定义为“Simpson”，即绝大多数情况下，我们在文献中看到的Simpson指数或者软件直接给出的Simpson指数结果，其实是Gini-Simpson指数，而并非经典Simpson指数。
四、均匀度（Evenness）
均匀度（Evenness）用于度量群落中相对物种丰富度。群落均匀度有多种表述方式，所代表的含义也略有区别，以Shannon均匀度和Simpson均匀度最为常见。
Shannon均匀度（Shannon’s evenness），又称Pielou均匀度（Pielou’s evenness），为群落实际的Shannon指数与具有相同物种丰富度的群落中能够获得的最大Shannon指数的比值。如果所有物种具有相同的相对丰度，则该值为1。
五、ACE指数（ACE）
ACE指数在生态学中同样作为度量物种丰富度的指标（Chao and Yang，1993）,其值越高代表群落物种越丰富。
在对某群落的研究中，我们获得了群落物种的观测数量，其中某些物种所被观测到的个体数量仅为1、2等很低的数值，通常表明这些物种在群落所有物种数量中的所占比例很低。但是如此低丰度的物种很容易由测量误差产生，从而出现较大的波动。为了排除这种测量误差所产生的干扰，我们指定一个稳定的阈值用于区分稀有物种和丰富物种，例如以10个个体为划分界限，将观测个体数小于或等于10的物种定义为稀有物种，大于10的即为丰富物种。与只观测到1、2等个体数量的物种相比，稀有物种的测量显得相对稳定。在群落物种测量中，根据丰富/稀有物种（根据指定阈值为界限）以及仅包含1个个体的物种，计算ACE指数用于估算群落中尚未被观测到的物种数量。其值越大，代表该群落中真实物种种类越多。
六、goods_coverage
goods_coverage指数的取值范围0~1，单位是一个概率。
goods_coverage常用在微生物16S/18S/ITS测序中，作为反映测序深度的指标。其值越接近于1，说明测序深度越合理，测序深度已经基本覆盖到样品中所有的物种；反之，测序深度不高，许多物种仅被测到了一次，暗示着很多低丰度物种可能尚未被测序测到。
七、系统发育多样性（Phylogenetic diversity）
谱系多样性（Phylogenetic diversity，PD）定义为在系统发育树上跨越给定的一组分类群所经过的所有系统发育分支的最小总长度（Faith，1992）。以下使用原作者文章中的一张图对谱系多样性作简要说明，对于更多详细内容，可参考原文章。
它同时考虑物种丰度以及进化距离，而上述几种类型的多样性指数是不关注进化距离的。
谱系多样性常用在微观群落的研究中。如基于16S/18S/ITS测序的微生物群落分析，在这些分析中进化树基于代表物种的16S/18S/ITS构建。

###################################################
###############R中计算alpha多样性###################
###################################################
`
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)`

`#定义函数
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')