Q:
什么是β多样性？
A:
β多样性是指在一个梯度上从一个生境到另一个生境所发生的多样性变化的速率和范围，它是研究群落之间的种多度关系。不同群落或某环境梯度上不同点之间的共有种越少，β多样性越大。精确地测定β多样性具有重要的意义。这是因为：①可以用来指示物种被生境隔离的程度；②可以用来度量生物多样性沿生境变化范围；③β多样性与α多样性一起构成了总体多样性或一定地段的生物异质性。
01
β多样性指数简介
β多样性指数可以分为群落成分指数与群落结构指数。
群落成分指数研究的是不同生境间物种数目的差异，例如最常用的Whittaker指数计算方法如下：
βw=S/mα-1
其中S为所研究系统中记录的物种总数；mα为各样方或样本的平均物种数。Whittaker指数越大表明不同生境的样方之间共有物种越少，物种多样性变化越大。对于微生物群落分析我们一般采用两两比较法，在Mothur中该指数的计算方法为：
其中ST为所有样品总物种数，SA和SB为样品A与B的样品数，A与B共有物种越少，Whittaker越大。而且由公式可以看出Whittaker值与总样品数目有关，若是三个样品两两比较，ST即为三个样品物种数目总和，所以所选样品数大于2时Mothur中计算的指数不一定是正数。
可以看出，成分指数仅考虑物种数，并没有考虑物种的相对丰度也即没有加权，然而微生物群落中微生物相对丰度差别很大，因此常用群落结构指数来分析。群落结构指数也叫生态距离，其计算方法有Euclidean、Manhattan、Bray Curtis、Jaccard等各种各样的计算方法，其中几种计算方法如下所示：
可以看出，欧氏距离即为n维空间2点之间直线距离，曼哈顿距离为各坐标值之差的绝对值之和，而切比雪夫距离为坐标值之差的绝对值的最大值。这些计算方法的缺点就是赋予不同物种相同的权重，也即无论是稀有物种还是优势物种相差1%的丰度距离相同，但是在生态学里由1%到2%和由91%到92%显然是不同的，因此在生态分析中群落数据常用的一种是Bray-Curtis指数，其计算方法如下所示：
也即两个样品之间的距离是每个物种丰度差值比上丰度之和，这时候显然由1%到2%距离要大于由91%到92%，但是有时候也会过分放大罕见物种的差别，可以去掉丰度过低的物种进行计算。
Mothur更多指数参见http://www.mothur.org/wiki/Calculators。
基于系统发育的β多样性一般通过OTU代表序列及其count table（也即OTU table）结合分析，使用unifrac方法计算生态距离，如下图所示：
对应的，其计算结果也分为考不虑物种丰度的UnweightedUnifrac以及考虑物种丰度的WeightedUnifrac。
02
β多样性距离计算
QIIME2、Mothur等软件均可计算β多样性距离指数。例如在Mothur中计算方法如下：
代码语言：javascript复制#计算样品两两之间的距离指数
summary.shared(shared=sample.0.03.0.03.subsample.shared, groups=QC1-QC3, calc=whittaker-braycurtis)
#计算样品间的距离矩阵
dist.shared(shared=sample.0.03.0.03.subsample.shared, calc=braycurtis, subsample=T, output=square)
#其中参数output=square则结果生成的是方形的矩阵，也即距离矩阵，可以通过设置output参数获得
#使用计算系统发育多样性产生的tre文件计算Unifrac距离矩阵
unifrac.unweighted(tree=sample.0.03.0.03.subsample.0.03.rep.phylip.tre, count=sample.0.03.0.03.subsample.0.03.rep.count_table, distance=square, random=F)
unifrac.weighted(tree=sample.0.03.0.03.subsample.0.03.rep.phylip.tre, count=sample.0.03.0.03.subsample.0.03.rep.count_table, distance=square, random=F)除了基于系统发育的Unifrac距离以外，微生物群落的距离矩阵均可以通过R计算获得。下面我们以生态学领域我们常用vegan包中的vegdist()函数为例，此函数使用方法如下所示：
代码语言：javascript复制vegdist(x, method="bray", binary=FALSE, diag=FALSE, upper=FALSE, ...)
#其中method的方法具体计算公式如下所示：
euclidean d[jk] = sqrt(sum(x[ij]-x[ik])^2),
binary：sqrt(A B-2*J),
where A and B are the numbers of species on compared sites, and J is the number of species that occur on both compared sites.
manhattan d[jk] = sum(abs(x[ij] - x[ik])),
binary：A B-2*J.
gower d[jk] = (1/M) sum(abs(x[ij]-x[ik])/(max(x[i])-min(x[i]))),
binary: (A B-2*J)/M,
where M is the number of columns (excluding missing values).
altGower d[jk] = (1/NZ) sum(abs(x[ij] - x[ik])),
binary: (A B-2*J)/(A B-J),
where NZ is the number of non-zero columns excluding double-zeros
canberra d[jk] = (1/NZ) sum ((x[ij]-x[ik])/(x[ij] x[ik])),
binary: (A B-2*J)/(A B-J),
where NZ is the number of non-zero entries.
bray d[jk] = (sum abs(x[ij]-x[ik]))/(sum (x[ij] x[ik])),
binary: (A B-2*J)/(A B).
kulczynski d[jk]=1-0.5*((sum min(x[ij],x[ik])/(sum x[ij]) (sum min(x[ij],x[ik])/(sum x[ik])),
binary: 1-(J/A J/B)/2.
morisita d[jk]=1-2*sum(x[ij]*x[ik])/((lambda[j] lambda[k]) * sum(x[ij])*sum(x[ik])),
binary: cannot be calculated,
where lambda[j] = sum(x[ij]*(x[ij]-1))/sum(x[ij])*sum(x[ij]-1).
Horn Like morisita, but lambda[j] = sum(x[ij]^2)/(sum(x[ij])^2),
binary: (A B-2*J)/(A B).
binomial d[jk] = sum(x[ij]*log(x[ij]/n[i]) x[ik]*log(x[ik]/n[i]) - n[i]*log(1/2))/n[i],
binary: log(2)*(A B-2*J),
where n[i] = x[ij] x[ik].
Cao d[jk] = (1/S) * sum(log(n[i]/2) - (x[ij]*log(x[ik]) x[ik]*log(x[ij]))/n[i]),
where S is the number of species in compared sites and n[i] = x[ij] x[ik].其中x为群落数据矩阵，其列名字为物种，行名字为样品；method为距离矩阵计算方法；binary为群落数据是否经过了有-无标准化；diag为是否显示对角线距离（对角线距离都是0）；upper为是否显示上三角部分，为TRUE的话计算结果为方阵。
最终距离的计算结果也要结合数据标准化处理（见1.4.2.1数据预处理）来进行评断，例如经过卡方转换后的数据使用欧氏距离方法计算会得到卡方距离矩阵。距离矩阵实际上代表的是对象之间的一种相异性（相似性），与数据标准化一样，距离矩阵只是一种数据转换方法，因此不需要进行假设检验。
我们可以基于PCoA比较相同群落不同距离计算对排序的影响，具体如下：
03
组间箱型图比较
对于一个样方内的样品点，或者一个处理组的样品，我们希望其群落相似也即距离相近，为此我们可以做组间或样方间β多样性箱线图，来比较观察每个样品组的样品之间的距离差异，具体如下所示：
代码语言：javascript复制#读取β多样性矩阵
beta=read.table("unweighted_unifrac_otu_table_13129.txt", header=TRUE, row.names=1)
#读取样品分组文件(样品名要与距离矩阵顺序一致)
group=read.table("groups.txt", header=TRUE, row.names=1)
group=factor(group[,1])
summary=summary(group)
#根据分组拆分距离矩阵
dist=as.matrix(beta)
group1=as.vector(dist[1:summary[1], 1:summary[1]])
group2=as.vector(dist[(summary[1] 1):(summary[1] summary[2]), (summary[1] 1):(summary[1] summary[2])])
group3=as.vector(dist[(summary[2] 1):(summary[2] summary[3]), (summary[2] 1):(summary[2] summary[3])])
group4=as.vector(dist[(summary[3] 1):(summary[3] summary[4]), (summary[3] 1):(summary[3] summary[4])])
weighted_unifrac=numeric(length=length(group1) length(group2) length(group3) length(group4))
weighted_unifrac[1:length(group1)]=group1
weighted_unifrac[(length(group1) 1):(length(group1) length(group2))]=group2
weighted_unifrac[(length(group2) 1):(length(group2) length(group3))]=group3
weighted_unifrac[(length(group3) 1):(length(group3) length(group4))]=group4
#去掉零值
weighted_unifrac[weighted_unifrac==0]=NA
groups=c(rep(names(summary)[1],length(group1)), rep(names(summary)[2],length(group2)), rep(names(summary)[3],length(group3)), rep(names(summary)[4],length(group4)))
data=data.frame(weighted_unifrac, groups)
mycol=c(99,81,503,562,76,96,495,52,619,453,71,134,448,548,655,574,36,544,89,120,131,596,147,576,58,429,386,122,87,404,466,124,463,552,147,45,30,54,84,256,100,652,31,610,477,150,50,588,621)
mycol=colors()[mycol]
par(mar=c(5,5,5,5))
boxplot(data$weighted_unifrac~data$groups, pch=" ", col=mycol, range=0.5, names=levels(data$groups), boxwex=0.5, ylab="Weighted_unifrac", xlab="Group", ylim=c(0.2, 0.9), cex.lab=1.2)作图结果如下所示：
每个样品组样品数目最好相同或相差不大，箱子越“扁”说明小组样品之间距离越近。
04
β多样性距离热图
样品间的β距离矩阵可以通过聚类距离热图来进行可视化，接下来我们均以weighted unifrac距离矩阵为例进行分析，方法如下所示：
代码语言：javascript复制dist=read.table("new.weighted.phylip.subsample.dist", header=FALSE)
rownames(dist)=dist[,1]
dist=dist[,-1]
colnames(dist)=t(rownames(dist))
dist=as.matrix(dist)
library(gplots)
otu_dist=as.matrix(dist)
mycol=colorRampPalette(c("yellow","orange", "red"))(50)
heatmap.2(dist, col=mycol, srtCol=45, trace="none", key.title=NA, dendrogram="both", cexRow=1, cexCol=1, keysize=1, key.ylab=NA, margins=c(5, 5), offsetRow=-0.3, offsetCol=-0.3, key.xlab="distance")作图结果如下所示：
—END—
