Q:
什么是群落测度？
A:
微生物群落的测度（measure）是指对群落矩阵数据的一种度量比较。测度可以用一系列指数（index）或系数（coefficient）来表示。对于单个对象（样品）的测度计算，可以采用α多样性指数来表示，而对于不同对象之间的比较，则可以采用β多样性指数或者距离。对于变量（物种或环境因子）之间的比较，则采用相关性来比较。群落测度的分析结果，可用于后续的排序分析、网络分析、聚类分析、判别分析等。
01
α多样性指数简介
α多样性指数反应群落内物种数量及其相对丰度，为群落内各物种利用同一生境互相竞争或共生的结果，比较不同样本的α多样性指数可以看出不同样本多样性差异。
α多样性指数可分为物种丰富度指数、物种均匀度指数、物种多样性指数三类：
①丰富度指数(communityrichness)反映的是群落内物种的丰富程度，常见的有sobs、chao、ace、jack、bootstrap等。sobs即根据测序结果在一定抽样水平实际观察到的物种数目（observedOTUs），chao、ace则分别用chao1和ace算法来估算一定抽样水平物种的总数目。
②均匀度指数(communityevenness)反映的是群落内物种的分配状况，也即各物种相对丰度的均匀程度，在同一群落内与丰富度指数呈正相关，常见的有simpsoneven、shannoneven、heip、smithwilson等。
③物种多样性指数(communitydiversity)反映的是物种丰富度和均匀度的综合状况，常见的有shannon-wienner（也即shannon）、simpson、invsimpson等。shannon指数反映的是物种丰度与均匀度，与这两者均呈正相关；simpson指数为在样本中抽取两条序列属于不同种的概率。
物种多样性指数为丰富度与均匀度的综合考察，其主要参数介绍如下。
① (Inv)Simpson指数
首先说一下反辛普森（Invsimpson）多样性指数计算公式如下：
其中ni为i物种（OTU）序列数目，N为总序列数，Sobs为观察到的物种数。公式的意思是在样本中随机抽取两条序列属于同一个物种（OTU）的概率，因此Invsimpson指数描述的是优势物种在群落中的作用和地位，也称为生态优势度，其值介于0和1之间。Invsimpson指数与其他多样性指数一般呈负相关，而Simpson指数描述的正好与Invsimpson指数相反，是随机抽取两条序列属于不同物种的概率，计算公式如下：
这两种指数之和为1，因此群落物种丰富度越高，每个物种相对丰度越均匀，则群落多样性越高，Simpson指数越高，Invsimpson指数越低。注意！不同软件计算的Simpson与Invsimpson指数其定义可能正好相反，在实际分析中需要具体问题具体分析。
② Shannon指数
香农（shannon）多样性指数计算公式如下：
从公式中可以看出，香农指数为正数。香农指数采用信息论的方法，反映的其实是群落的不确定性。群落多样性越高，物种越丰富、均匀度越高，则香农指数越高。更多多样性指数请参见http://www.mothur.org/wiki/Calculators。
02
α多样性指数计算
扩增子数据分析常用的工具QIIME2和Mothur均能计算α多样性指数。在Mothur中summary.single指令进行计算分析，产生多个样品的各种多样性指数计算结果，计算方法如下：
代码语言：javascript复制summary.single(shared=sample.0.03.0.03.subsample.shared, calc=sobs-chao-shannon-simpson, subsample=T)因为序列总数不一样，测序深度就会有所不同，因此要使得不同样品的多样性指数计算有意义，指数计算必须在统一抽样水平，我们可以设置参数subsample为具体的序列数目，设置subsample=T则Mothur会自动将序列数最少的样品的序列数作为统一的抽样水平（若是已经进行过抽平可以不设置这个参数），运行结束后会产生summary文件，展示每个样品各种指数计算结果，如下图所示:
其中chao为平均值，chao_lci为最低值，chao_hci为最高值，以此类推。
基于OTU的代表序列，我们开可以计算系统发育多样性指数，方法为计算样品中系统发育树总分支长度。首先计算代表序列的距离矩阵并建树：方法如下所示：
代码语言：javascript复制dist.seqs(fasta=sample.0.03.0.03.subsample.0.03.rep.fasta, output=lt, processors=20)
clearcut(phylip=sample.0.03.0.03.subsample.0.03.rep.phylip.dist)注意所输入的fasta文件也必须是比对的文件。然后进行多样性指数计算：
代码语言：javascript复制phylo.diversity(tree=sample.0.03.0.03.subsample.0.03.rep.phylip.tre, count=sample.0.03.0.03.subsample.0.03.rep.count_table)运行结束后，会生成summary文件，里面为每个样品的系统发育多样性指数，如下所示：
在R中可以利用vegan包中的diversity()函数来计算微生物群落的alpha多样性指数，但不能计算系统发育多样性指数。
03
组间箱型图比较
在大规模环境微生物测序中，样品间多样性比较太冗杂，我们一般根据实际情况对样品进行分组，将Mothur计算的多样性指数分组统计进行比较，我么可以采用箱形图来进行展示。如下图所示一共有30个样品的香农指数数据：
如果直接做图分析，规律性很差，现根据pH分为A、B、C三个组，做一组箱型图来比较不同多样性指数的变化规律，方法如下所示：
代码语言：javascript复制alpha=read.table("diversity.txt",head=T)
rownames(alpha)=alpha[,1]
alpha=alpha[,-1]
par(mfrow=c(2,2))
boxplot(alpha$sobs~alpha$Group_ID, pch=" ", col="lightblue", range=0.8, names=c("pH>5", "5>pH>3", "pH<3"), boxwex=0.5, ylab="Observed OTUs")
boxplot(alpha$shannon~alpha$Group_ID, pch=" ", col="lightblue", range=0.8, names=c("pH>5", "5>pH>3", "pH<3"), boxwex=0.5, ylab="Shannon Index")
boxplot(alpha$simpson~alpha$Group_ID, pch=" ", col="lightblue", range=0.8, names=c("pH>5", "5>pH>3", "pH<3"), boxwex=0.5, ylab="Simpson Index")
boxplot(alpha$phyloDiversity~alpha$Group_ID, pch=" ", col="lightblue", range=0.8, names=c("pH>5", "5>pH>3", "pH<3"), boxwex=0.5, ylab="Phylogenetic Diversity")其中range=1设置了触须伸出箱子的最大长度，range=0则触须为最大与最小值。触须外的离群点用“ ”显示出来，这时候我们可以看出pH<3的强酸性条件下微生物多样性大大降低。作图结果如下所示：
04
α多样性稀释曲线
稀释曲线展示的是在不同测序深度（即抽样水平）下群落多样性的变化。多样性指数的稀释曲线可以使用Mothur中的rarefaction.single指令进行分析产生，rarefaction.single为采用不重复抽样的方法产生多样性指数随序列抽样水平变化的稀释曲线，这里的多样性指数一般为丰度指数，也可以是shannon等多样性指数；而phylo.diversity命令则可以产生系统发育多样性的稀释曲线。方法如下所示：
代码语言：javascript复制rarefaction.single(shared=sample.0.03.0.03.subsample.shared, calc=sobs-chao-shannon-simpson, freq=100, processors=20)
phylo.diversity(tree=sample.0.03.0.03.subsample.0.03.rep.phylip.tre, count=sample.0.03.0.03.subsample.0.03.rep.count_table, rarefy=T, freq=10)其中freq为抽样频率参数，freq=100也即抽样间隔为100条序列（1、100、200、300）。运行完成后每个指数均生成一个结果文件，其中sobs结果为rarefaction文件。打开生成的rarefaction文件如下图所示，每个样品的多样性指数第一列数据为平均值，lci为最低值，hci为最高值：
我们可以简单处理下这些结果文件，然后使用R进行作图，方法如下所示：
代码语言：javascript复制alpha=read.table("new.subsample.rarefaction", header=TRUE)
mycol=c(119,132,147,454,89,404,123,463,552,28,54,84,256,100,558,43,652,31,610,477,588,99,81,503,562,76,96,495)
mycol=colors()[mycol]
n=max(alpha[, 1])
m=max(alpha[, -1])
layout(1, 1, 1)
par(mar=c(5,5,5,5))
plot(alpha[, 1], alpha[, 2], type="l", xlab="Sequences sampled", ylab="Observed OTUs", xlim=c(0, 1.1*n), ylim=c(0, 1.1*m), col=mycol[1] , tcl=0.2, lwd=2, font.lab=2)
for (i in 3:ncol(alpha)) {
points(alpha[, 1], alpha[, i], type="l", col=mycol[i-1], lwd=2)
}
for (i in 2:ncol(alpha)) {
x=alpha[, 1]
y=alpha[, i]
text(x[nrow(alpha)] 500, y[nrow(alpha)], labels=colnames(alpha)[i], cex=1, col=mycol[i-1])
}作图结果如下所示：
稀释曲线可以反映测序深度，如果曲线末端与X轴接近平行，说明测序数据足够反映真实的群落信息，曲线最后的指数值也是可信的。虽然sobs的稀释曲线主要是针对丰富度指数进行计算，然而稀释曲线更多的反映多样性信息。假如两个样本丰富度指数相同，但其均匀度不同，则其中间抽样过程中产生的丰富度指数是不同的，也即其稀释曲线虽然两端重合，但形状不同。
—END—
