0 导语

时间过得真快。
写完《群落多样性之Beta多样性（一）》还是2019年高考前夕，虽说今年高考拖延一个月，但也早已经结束，当年的考生们已经大二了。
严重的拖延症！！！
回顾一下，（一）中比较粗暴地讲了Beta多样性的两个重要部分：
1）群落间多样性距离的计算；
2）距离的展示，降维思想，即PCA、PCoA和NMDS等。
简言之，就是距离的计算和展示。
Beta多样性（二）打算深入细节对距离的计算展开介绍，如果大哥们只喜欢粗暴去理解原理，不感兴趣细节的话，只须阅读（一）即可。

1 相似度和距离

相似度（Similarity）和距离（Distance）是一对反义的概念，一般来说，相似度越高，距离越近。两者在模式分析问题中起着关键作用，如分类、聚类等。
相似度/距离计算的输入数据类型一般分为两种：布尔型（即真和假）和数值型，其对应的距离分别叫做非加权距离和加权距离。
非加权距离是依据特征属性在样本中体现的真假来计算，真即1，假即0，这属于定性资料，并不考虑特征属性的定量信息。放在宏基因组分析中，特征属性即为物种，即某物种是否存在于样本中，不存在即为0，存在即为1。
数值型数据是特征属性的在样本中的定量信息，属于定量资料，计算过程中考虑特征属性的定量信息。在宏基因组分析中，加权距离的计算不仅考虑物种的存在与否还会考虑物种的丰度。

对于非加权距离的计算，我们需要知道物种在样本 $i$ 和 $j$ 中存在与否的各个集合信息（表1）。

表1 物种在双样本（ $i$ 和 $j$ ）中是否存在的各个集合汇总表

j\i	1 (Presence)	0 (Absence)	Sum
1 (Presence)	$a=i*j$	$b=\bar{i}*j$	$a+b$
0 (Absence)	$c=i*\bar{j}$	$d=\bar{i}*\bar{j}$	$c+d$
Sum	$a+c$	$b+d$	$n=a+b+c+d$

如表1中表示， $a$ 表示在样本 $i$ 和 $j$ 中都存在即都为真即 $(1,1)$ 的物种数量； $b$ 表示在样本 $i$ 和 $j$ 中分别为一假一真 $(0,1)$ 的物种数量；以此类推， $c$ 则表示布尔值为 $(1,0)$ 的物种数量； $d$ 则表示 $(0,0)$ 的物种数量。 $n$ 为物种的总数量。

非加权的距离计算，理解这张表很重要，因为马上在下文就能体现出其重要性。
在实际的应用中，不同类型数据性能的评估往往更加依赖一个更加独特且合适的相似度/距离计算方式，因此各领域的科学家和工程师们花了很多精力在来寻找最有意义的相似性/距离度量方法。图1 展示了各种距离算法的发展时间轴。

图1 二元相似性/距离度量方法发展年表[1]

Choi等[1] 总结了2005年之前所有已发表的非加权距离计算公式（基于表1），共有76个之多（表2）。表2由于篇幅限制展示了部分非加权距离计算公式，如对这些计算方法感兴趣，请阅读Choi等发表的文章的原文。

表2 基于二进制数据的距离计算公式[1]

公式	编号
$S_{JACCARD}=\frac{ a} {a+b+c}$	$(1)$
$S_{DICE}=\frac{2a}{2a+b+c}$	$(2)$
$S_{CZEKANOWSKI}=\frac{2a}{2a+b+c}$	$(3)$
$S_{3W-JACCARD}=\frac{3a}{3a+b+c}$	$(4)$
$S_{NEI\&LI}=\frac{2a}{(a+b)+( a+c)}$	$(5)$
$......$	$......$
$S_{BARONI-URBANI\&BUSER-I}=\frac{\sqrt{ad}+a}{\sqrt{ad}+a+b+c}$	(71)
$S_{BARONI-URBANI\&BUSER-II}=\frac{\sqrt{ad}+a-(b+c)}{\sqrt{ad}+a+b+c}$	(72)
$S_{PEIRCE}=\frac{ab+bc}{ab+2bc+cd}$	(73)
$S_{EYRAUD}=\frac{n^2(na-(a+b)(a+c))}{(a+b)(a+c)(b+d)(c+d)}$	(74)
$S_{TARANTULA}=\frac{\frac{a}{a+b}}{\frac{c}{c+d}} =\frac{a(c+d)}{c(a+b)}$	(75)
$S_{AMPLE}=\mid{}\frac{\frac{a}{a+b}}{\frac{c}{c+d}}\mid{}=\mid{}\frac{a(c+d)}{c(a+b)}\mid{}$	(76)

由于计算距离的公式众多，在此不能一一展示。
本文着重对宏基因组中应用较为广泛的Bray&Curtis和Unifrac距离展开详细的解读。
上述丰度计算方法尽管是只考虑OTU在样本中存在与否的二元相似度和距离的计算公式，但是通常情况下稍做转换即可得出基于OTU定量信息的计算方法。

2 Bray&Curtis距离的计算方法

Bray&Curtis距离以开发者J. Roger Bray和John T. Curtis的名字命名[2]，可基于物种的有无和物种的丰度两种方式比较两个群落微生物的组成差异。尽管我还没发现有这么叫的，但我这里依然分别称之为非加权和加权Bray&Curtis距离。

非加权Bray&Curtis距离计算公式如下[1]：
$D_{Bray\&Curtis}=1-\frac{2a}{(a+b)+(a+c)}=\frac{b+c} {2a+b+c}$
上式中， $1$ 表示标准化的全集合， $2a$ 表示 $i$ 样本和 $j$ 样本中共有的OTU数量的2倍，由于这类OTU需要在每个样本都要统计一次因此需要乘以2， $(a+b)+(a+c)$ 表示 $i$ 样本和 $j$ 样本中的OTU数量之和。整个分式 $\frac {2a} {(a+b)+(a+c)}$ 则是表示共有的OTU占总OTU中的比例，其意义则是共有OTU丰度占总OTU丰度的比例，即相似度的一种度量，而距离=1-相似度。
比如有两个样本A和B，OTU及其对应丰度信息如表3所示。

表3 样本A和B的OTU丰度表

Sample	OTU1	OTU2	OTU3	OTU4	OTU5
A	10	2	0	5	0
B	0	8	3	5	2

此时不考虑丰度，参考表1分别计算各个集合的大小， $a=2$ ， $b=2$ ， $c=1$ ， $d=0$ 。
代入公式：
$D_{Bray\&Curtis}=\frac{b+c}{2a+b+c}=\frac{2+1}{2\times2+1+1}=0.5$

加权Bray&Curtis距离计算公式如下[3]：

$D_{Bray\&Curtis}=1-2\frac{\sum min(S_{A,i},S_{B,i})}{\sum S_{A,i}+\sum S_{B,i}}$

上式中， $S_{A,i}$ 和 $S_{B,i}$ 分别表示第 $i$ 个OTU在样本A或样本B中的丰度。 $min(S_{A,i}, S_{B,i})$ 为取第 $i$ 个OTU在两个样本中丰度的最小值（表1），分子 $\sum min(S_{A,i}, S_{B,i})$ 本质上是对存在于两个样本中OTU的进行求和，乘以2的目的是因为两个样本的共有的OTU在每个样本中各自计算1次，分母 $\sum S_{A,i}+\sum S_{B,i}$ 则是两个样本中所有OTU个数之和。

2.1和2.2中的两个公式本质上的意义是一致的，唯一的不同就是第一个公式是基于OTU的存在与否，而第二个公式是基于OTU的丰度来计算两个群落微生物的组成差异。

3 UniFrac距离

Bray&Curtis等距离的计算方法往往基于一个前提，就是假设属性之间是独立的，而事实上往往这些属性是存在一定关系的，因此计算距离的时候也可以把这些属性之间的关系考虑在内。在宏基因组的研究中，研究者可以选择用UniFrac距离来兼顾物种之间的进化关系计算群落间的距离。
UniFrac距离算法，诞生于2005年底，由Lozupone和Knight开发[4]。
我们可以顾名思义地猜测一下两位大神作者的取UniFrac这个名字的意图。
Uni，即Unique，独有的；Frac，即Fraction，分数，也可以看成是比例。
那么具体是如何Uni，又是怎么个Frac呢？
请诸位大哥继续往下看。
UniFrac距离包括非加权（Unweighted）和加权（Weighted）UniFrac距离，这里只介绍非加权UniFrac距离的原理和计算方法。

3.1 非加权UniFrac距离

由于UniFrac距离考虑物种进化关系，因此在计算UniFrac距离之前，我们要先计算其OTU之间的系统发育关系。
如图2所示，已知两个微生物群落样本，分别由红色方块和绿色圆形表示。紫色的分支表示来自某单一群落的后代。UniFrac用紫色树枝的长度除以树的总树枝长度来衡量这两个群落之间的进化差异。左边的例子显示了两个群落的进化差异很小，因此UniFrac距离相对较低。相反，右边的例子显示了两个最大差异的群落，因为它们最低的共同祖先是树的根[3]。

图2 Unweighted UniFrac距离计算方法图解

非加权UniFrac距离可使用下面的公式计算[5]。
$Unweighted\_UniFrac=\frac{\sum_{i=1}^{N}l_{i}|A_{i}-B_{i}|}{\sum_{i=1}^{N}l_{i}max(A_{i},B_{i})}$

其中， $N$ 是在树上的节点数量， $l_i$ 是节点 $i$ 及其父节点之间的长度， $A_i$ 和 $B_i$ 取值 $0$ 或 $1$ 分别表示节点 $i$ 的子节点是否存在群落 $A$ 和群落 $B$ 。

UniFrac可以应用于多个群落的距离计算，具体方法是：先构建一个由所有序列组成的树，任何一对群落之间的距离则通过修剪掉其他群落的树枝计算，这样树就只包含来自这对群落的序列。（这在上面给出的公式中是隐式的，因为没有出现在任何一对群落中的节点的统计量为零，所以不影响结果。）

接下来应用随机排列检验计算两个群落之间的距离是否显著大于随机预期。如图3所示，这是通过将群落标签分配到序列中，同时保持树不变来实现的。假设检验获得的 $p$ 值是导致UniFrac距离大于或等于观测数据的随机排列的概率（图3）。

图3 随机排列检验检验群落距离显著性

两个群落的距离是这么计算的，那么三个的呢？
上述仅介绍了两个群落之间的距离计算方法，图4较为清晰地展示了三个群落计算非加权UniFrac的具体操作流程，即A和B构建系统发育树计算UniFrac距离，C为随机排列检验，D则是距离矩阵的计算。

图4 计算Unweighted UniFrac距离矩阵的基本流程。正方形、三角形和圆形表示来自不同群落的序列。节点分支如果全部处于单一群落，则为黑色；如果是节点分支是多群落共享共享的，则为灰色。

3.2 加权UniFrac距离的计算

加权UniFrac算法实际上可以看作是非加权UniFrac的扩展版本，进一步引入了属性的量化信息。
下图展示了将加权UniFrac应用于两个由红色方块和绿色圆圈表示的的两个群落样本的计算结果。对各个分枝进行加权，并计算各群落中分枝下的序列的相对比例。紫色的分支是那些导致分子非零值的分支。考虑到序列之间的演化速率差异，加权UniFrac度量被标准化，使得每个序列对计算的距离贡献相等。

图5 加权UniFrac距离的计算

上图的公式中， $N$ 是在树上的节点数量， $S$ 是树所表示的序列的数目， $l_i$ 是节点 $i$ 及其父节点之间的长度， $Lj$ 是从树根到树梢即序列 $j$ 的总分支长度， $A_i$ 和 $B_i$ 是分别来自群落 $A$ 和 $B$ 的节点数量，而 $A_T$ 和 $B_T$ 分别是群落样本 $A$ 和 $B$ 序列的总数。

加权的UniFrac可以应用于多个群落样本距离矩阵的计算，构建一个由所有序列组成的树。每对群落样本之间的距离是通过修剪树来计算的，这样树就只包含来自这对群落的序列，可以通过随机排列测试来测试两个社区之间的距离是否大于预期，这与未加权的UniFrac类似。

3.2 加权和非加权的抉择

该算法作者在2006年的一篇文章中[6]告诉我们对数据集同时进行加权UniFrac和非加权UniFrac算法往往会导致不同的结果，因为它们会导致关于微生物群落相似性和驱动微生物多样性的因素的不同结论。
一般来说，如果群落结构物种的种类差别较大的情况下，用非加权UniFrac距离要比用非加权UniFrac距离更能找到区别；反之如果仅仅是物种丰度差别大的话，则选择加权Unifrac距离。

3.3 Fast UniFrac算法

2010年，Hamady等[7]在Unifrac算法的基础上开发了 Fast UniFrac算法。如图4所示，以四个群落样本（R、B、O和G）为例，
1）首先，将各个群落样本进行根据每个OTU的有无进行编码，比如0号OTU仅仅存在于R中，则编码为1000；1号OTU仅存在于B中，则编码为0100，……以此类推。
2）接下来，去掉暂时不相关的群落，选取要比较的群落和对应的编码。
3）使用并集算法对环境集进行比较，将状态分配给每个内部节点，比如5号和6号OTU分别编码为00和01，那么它们共同的上一级父节点11号节点则编码为01；
4）最后，计算各个群落样本的枝长百分比再求和，就是非加权Unifrac距离，实际上图4h中的公式与3.1 非加权Unifrac距离是等效的。

图4 UniFrac和Fast UniFrac在流程上的差异（非加权）。UniFrac算法：a）环境以集合的形式存储在树对象中；b）这棵树被修剪，只包括那些通向想要的环境的枝干；c）使用集合算法对环境集进行比较，将状态分配给每个内部节点；d）结果由另一个树遍历计算。Fast UniFrac算法：e）环境被存储为提示x环境计数的数组；f）通过对该数组进行切片来选择所选环境；g）使用数组片上的数组操作计算内部状态；h）将关联数组和连接两个环境中的一个或两个环境的节点的分支长度的乘积相加，从而计算UniFrac值。基于数组的方法可以大大提高效率。

4 写在后面

最后，想说一说一个文学作品中永恒的主题——爱情。
爱情在文学作品中，可以说是突破性非常强的一个东西。
爱情可以突破社会阶层，突破地域，突破宗教信仰，突破年龄……
作为生物狗，到这里我都还可以理解。
不过也有我无法理解，因为他们跨越了物种。
很明显，这不科学。

若干年前，有一首歌火遍大江南北，一首让我乍一听就满脸鸡皮疙瘩的歌。我一说名字你准知道——《狼爱上羊》，内容通篇如歌名，没听过的可以搜一下。

你要说马和驴之间的爱情，我还可以勉强接受。
虽然染色体数量不一样，无法产生可育后代即骡子，但他们毕竟长相差别不大，也难保马或驴不脸盲啊。
狼和羊本是食物链关系，他们是不会有结果的。
不过，还有更过分的。

曾经，我的某室友QQ昵称叫“飞鸟和鱼”。
啥意思？这令我很费解，如果食物链，那么为啥不叫“水鸟和鱼”？
后来才知道，原来是有这么一首诗《飞鸟和鱼》，诗里面有这么一句：
“世界上最遥远的距离，不是我站在你面前，你却不知道我爱你，
而是明明知道彼此相爱，却不能在一起。”通篇都是这类句型。

又是一对食物链的爱情，从生物分类学角度来讲，狼和羊尚同属哺乳纲，飞鸟和鱼分属鱼纲和鸟纲，已经跨纲了……然而，跨纲算个屁，还有跨界的。你听（此处应有歌声）：

“……我爱你爱着你就像老鼠爱大米……”
这歌一定听过吧，是不是很熟悉？
又是食物链的爱情……。

人类社会中尚且不允许兄妹之间产生爱情，何况不同物种？
如此，可以总结一下，突破了物种的爱情是不会有结果的。
我们且看演化生物学家恩斯特·麦尔的定义：“物种是能够相互配育的自然种群的类群，这些类群与其它这样的类群在生殖上相互隔离着。”

这些距离，由远及近，总结起来很有层次感。
最远的，群落层面，Beta多样性；
近一点呢，物种间层面，物种和物种之间的距离，不同物种之间性交不可育；
再近一点，物种内层面，即物种间各个个体之间的距离，物种内性交可育。
还有番外篇，个体内，基因与基因间的按重组率遗传距离，以厘摩（cM）为单位，做过遗传图谱的大哥们都知道。
以上这些均可通过DNA序列相关数据计算获得。

一般来说，两个对象之间的距离很直观，但是意义极为有限，甚至可以说没有意义。为了让这个距离显得有意义，现实的分析过程中，往往需要针对多个对象，分别计算每两个对象之间的距离，最后得到一个对称的距离矩阵（如表4）。尽管应用该矩阵可以得知两个对象之间的距离，但是非常不直观，因而难以对这些对象进行分类。至于如何更直观地对距离进行可视化展示，我们就留到后面再说吧。
表4 距离矩阵

Object	A	B	C	D	E
A	0	2	7	5	8
B	2	0	3	5	2
C	7	3	0	2	3
D	5	5	2	0	6
E	8	2	3	6	0

备注：本文于2020年11月22日发表于e媛微生态公众号。

参考文献
[1] Choi S-S, Cha S-H, Tappert C C: A survey of binary similarity and distance measures[J]. Journal of Systemics, Cybernetics and Informatics,2010, 8(1):43-48.
[2] Bray J R, Curtis J T: An ordination of the upland forest communities of southern Wisconsin[J]. Ecological monographs,1957, 27(4):326-349.
[3] https://mothur.org/wiki/braycurtis/
[4] C L, R K: UniFrac: a new phylogenetic method for comparing microbial communities.[J]. Appl Environ Microbiol,2005, 71(12):8228.
[5] https://mothur.org/wiki/unweighted_unifrac_algorithm/
[6] Lozupone C A, Hamady M, Kelley S T, et al: Quantitative and qualitative β diversity measures lead to different insights into factors that structure microbial communities[J]. Appl Environ Microbiol, 2007, 73(5): 1576-1585.
[7] https://mothur.org/wiki/weighted_unifrac_algorithm/
[8] Hamady M, Lozupone C, Knight R: Fast UniFrac: facilitating high-throughput phylogenetic analyses of microbial communities including analysis of pyrosequencing and PhyloChip data[J]. ISME J, 2010, 4(1):17-27.

群落多样性之Beta多样性（二）