物种数量及多样性的外推
物种数量及多样性的外推
前言
对于微生物群落的分析，最基本的就是得到群落的物种数量，也即OTU的数量。在此基础上，通过观察到的物种数量进行合理的外推，可以获得理论物种数。对于这两种物种数量的计算，R中都可以非常方便的完成。本文对此进行介绍。
计算观测到的物种
先构建一个OTU表：
>otu = matrix(sample(1:100, 25, replace = TRUE), nrow = 5, ncol = 5)
>rownames(otu) <- paste0("OTU", 1:nrow(otu))
>colnames(otu) <- paste0("Sample", 1:ncol(otu))
>otu
Sample1 Sample2 Sample3 Sample4 Sample5
OTU1 56 84 97 56 24
OTU2 37 92 48 95 86
OTU3 84 57 68 74 15
OTU4 57 8 66 53 91
OTU5 40 36 44 66 65计算物种数，可以通过求每一列大于零的个数直接得到，即
>S = colSums(otu>0);S
Sample1 Sample2 Sample3 Sample4 Sample5
5 5 5 5 5 也可以由vegan包中的函数求得：
>library(vegan)
>S <- specnumber(t(otu));S
Sample1 Sample2 Sample3 Sample4 Sample5
5 5 5 5 5 对于物种数量的外推，目前使用较为广泛的有Chao1,Chao2, ACE, ICE, Jack1, Bootstrap等指数。
Chao1, ACE是基于物种abundance信息对物种数进行的外推； Chao2, ICE是基于物种incidence信息对物种数进行的外推。会先把abundance数据转化为0-1数据再进行计算。 Jack1和Bootstrap方法根据其参数不同分别可实现基于abundance或incidence的外推。
其公式可以参考（阅读原文可查看）：
http://viceroy.eeb.uconn.edu/estimates/estimatespages/estsusersguide/estimatesusersguide.htm
简单的来说，Chao的方法主要基于singletons和doubletons来计算。ACE和ICE主要基于singletons和出现次数小于等于10次的稀有物种来计算。在对群落进行抽样的时候，如果还存在没有被发现的新物种，那么就一直会观察到低丰度稀有物种的出现。因此虽然每种方法的计算公式不同，但是核心思想都是利用稀有物种的个数来对群落整体物种数进行估计。因此，生成OTU时的方法也会显著影响对物种估计的值。
计算上我常用两个包：fossil和iNEXT。
fossil: Palaeoecological and Palaeogeographical Analysis Tools
其计算上述多样性指数的函数也非常简单，如：
>library(fossil)
>chao1(otu)
>chao2(otu)
>ACE(otu)
>ICE(otu)
>jack1(otu) #基于abundance
>jack1(otu,abund = FALSE) #基于incidence
>bootstrap(otu,abund = FALSE) #基于abundance
>bootstrap(otu, abund = TRUE,) #基于incidence
#spp.est 估计物种多样性
#spp.est(x, rand = 10, abund = TRUE, counter = FALSE, max.est = 'all')
#rand计算次数。counter计算过程是否显示。max.est参与计算的样本数
#abundance数据：计算richness,chao1,ACE,jack1及95%上下区间
#incidince数据：计算richness,chao2,ICE,jack1及95%上下区间
>spp.est(otu, abund = TRUE, counter = TRUE)
N.obs S.obs S.obs(+95%) S.obs(-95%) Chao1 Chao1(upper) Chao1(lower) ACE ACE(upper)
[1,] 1 5 5 5 5 5 5 5 5
[2,] 2 5 5 5 5 5 5 5 5
[3,] 3 5 5 5 5 5 5 5 5
[4,] 4 5 5 5 5 5 5 5 5
[5,] 5 5 5 5 5 5 5 5 5
ACE(lower) Jack1 Jack1(upper) Jack1(lower)
[1,] 5 5 5 5
[2,] 5 5 5 5
[3,] 5 5 5 5
[4,] 5 5 5 5
[5,] 5 5 5 5iNEXT: iNterpolation and EXTrapolation of Hill number
这个包就是Anne Chao等人开发的，主要用于计算Hill number。 关于Hill number之前有文章介绍过，详见：
拓展种-面积关系(SAR)为多样性-面积关系(DAR)
对于物种数量的计算同时也具有abundance和incidence两种方法，另外该包还可以计算shannon，simpson及他们的外推估计值。
#计算物种数和估计的物种数
>library(iNEXT)
>ChaoRichness(spider$Girdled, datatype="abundance")
Observed Estimator Est_s.e. 95% Lower 95% Upper
Sample1 5 5 0.000 5 5.000
Sample2 5 5 0.018 5 5.036
Sample3 5 5 0.000 5 5.000
Sample4 5 5 0.000 5 5.000
Sample5 5 5 0.001 5 5.001
#计算shannon和估计的shannon
>ChaoShannon(spider$Girdled, datatype="abundance")
Observed Estimator Est_s.e 95% Lower 95% Upper
Sample1 1.565 1.572 0.017 1.565 1.605
Sample2 1.421 1.428 0.030 1.421 1.486
Sample3 1.569 1.575 0.016 1.569 1.605
Sample4 1.586 1.592 0.012 1.586 1.615
Sample5 1.433 1.440 0.031 1.433 1.501
#simpson
>ChaoSimpson(spider$Girdled, datatype="abundance")
Observed Estimator Est_s.e. 95% Lower 95% Upper
Sample1 0.781 0.784 0.008 0.781 0.800
Sample2 0.738 0.740 0.010 0.738 0.760
Sample3 0.783 0.786 0.007 0.783 0.800
Sample4 0.790 0.793 0.005 0.790 0.803
Sample5 0.738 0.740 0.010 0.738 0.760
##Hill的计算和外推
>data(spider)
>out1 <- iNEXT(spider, q=0, datatype="abundance") #q为hill的阶数，可选0，1，2
>out1$AsyEst # species,shannon及simpson的观测值和预测值。这个是主要需要的数据。
Site Diversity Observed Estimator s.e. LCL UCL
1 A Species richness 5.000 5.000 0.000 5.000 5.000
2 A Shannon diversity 4.782 4.817 0.082 4.782 4.978
3 A Simpson diversity 4.575 4.636 0.158 4.575 4.946
4 B Species richness 5.000 5.000 0.018 5.000 5.036
5 B Shannon diversity 4.140 4.171 0.133 4.140 4.431
—END—