0 導(dǎo)語

時間過得真快着逐。
寫完《群落多樣性之Beta多樣性（一）》還是2019年高考前夕，雖說今年高考拖延一個月同诫，但也早已經(jīng)結(jié)束，當(dāng)年的考生們已經(jīng)大二了。
嚴(yán)重的拖延癥！Ｅで恪！
回顧一下挽绩，（一）中比較粗暴地講了Beta多樣性的兩個重要部分：
1）群落間多樣性距離的計算膛壹；
2）距離的展示，降維思想唉堪，即PCA恢筝、PCoA和NMDS等。
簡言之巨坊，就是距離的計算和展示。
Beta多樣性（二）打算深入細(xì)節(jié)對距離的計算展開介紹此改，如果大哥們只喜歡粗暴去理解原理趾撵，不感興趣細(xì)節(jié)的話，只須閱讀（一）即可共啃。

1 相似度和距離

相似度（Similarity）和距離（Distance）是一對反義的概念占调，一般來說，相似度越高移剪，距離越近究珊。兩者在模式分析問題中起著關(guān)鍵作用，如分類纵苛、聚類等剿涮。
相似度/距離計算的輸入數(shù)據(jù)類型一般分為兩種：布爾型（即真和假）和數(shù)值型言津，其對應(yīng)的距離分別叫做非加權(quán)距離和加權(quán)距離。
非加權(quán)距離是依據(jù)特征屬性在樣本中體現(xiàn)的真假來計算取试，真即1悬槽，假即0，這屬于定性資料瞬浓，并不考慮特征屬性的定量信息初婆。放在宏基因組分析中，特征屬性即為物種猿棉，即某物種是否存在于樣本中磅叛，不存在即為0，存在即為1萨赁。
數(shù)值型數(shù)據(jù)是特征屬性的在樣本中的定量信息弊琴，屬于定量資料，計算過程中考慮特征屬性的定量信息位迂。在宏基因組分析中访雪，加權(quán)距離的計算不僅考慮物種的存在與否還會考慮物種的豐度。

對于非加權(quán)距離的計算掂林，我們需要知道物種在樣本 $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)$ 的物種數(shù)量泻帮； $b$ 表示在樣本 $i$ 和 $j$ 中分別為一假一真 $(0,1)$ 的物種數(shù)量精置；以此類推， $c$ 則表示布爾值為 $(1,0)$ 的物種數(shù)量锣杂； $d$ 則表示 $(0,0)$ 的物種數(shù)量脂倦。 $n$ 為物種的總數(shù)量。

非加權(quán)的距離計算元莫，理解這張表很重要赖阻，因為馬上在下文就能體現(xiàn)出其重要性。
在實際的應(yīng)用中踱蠢，不同類型數(shù)據(jù)性能的評估往往更加依賴一個更加獨特且合適的相似度/距離計算方式火欧，因此各領(lǐng)域的科學(xué)家和工程師們花了很多精力在來尋找最有意義的相似性/距離度量方法。圖1 展示了各種距離算法的發(fā)展時間軸茎截。

圖1 二元相似性/距離度量方法發(fā)展年表[1]

Choi等[1] 總結(jié)了2005年之前所有已發(fā)表的非加權(quán)距離計算公式（基于表1）苇侵，共有76個之多（表2）。表2由于篇幅限制展示了部分非加權(quán)距離計算公式企锌，如對這些計算方法感興趣榆浓，請閱讀Choi等發(fā)表的文章的原文。

表2 基于二進制數(shù)據(jù)的距離計算公式[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)

由于計算距離的公式眾多撕攒，在此不能一一展示陡鹃。
本文著重對宏基因組中應(yīng)用較為廣泛的Bray&Curtis和Unifrac距離展開詳細(xì)的解讀烘浦。
上述豐度計算方法盡管是只考慮OTU在樣本中存在與否的二元相似度和距離的計算公式，但是通常情況下稍做轉(zhuǎn)換即可得出基于OTU定量信息的計算方法杉适。

2 Bray&Curtis距離的計算方法

Bray&Curtis距離以開發(fā)者J. Roger Bray和John T. Curtis的名字命名[2]谎倔，可基于物種的有無和物種的豐度兩種方式比較兩個群落微生物的組成差異。盡管我還沒發(fā)現(xiàn)有這么叫的猿推，但我這里依然分別稱之為非加權(quán)和加權(quán)Bray&Curtis距離片习。

非加權(quán)Bray&Curtis距離計算公式如下[1]：
$D_{Bray\&Curtis}=1-\frac{2a}{(a+b)+(a+c)}=\frac{b+c} {2a+b+c}$
上式中， $1$ 表示標(biāo)準(zhǔn)化的全集合蹬叭， $2a$ 表示 $i$ 樣本和 $j$ 樣本中共有的OTU數(shù)量的2倍藕咏，由于這類OTU需要在每個樣本都要統(tǒng)計一次因此需要乘以2， $(a+b)+(a+c)$ 表示 $i$ 樣本和 $j$ 樣本中的OTU數(shù)量之和秽五。整個分式 $\frac {2a} {(a+b)+(a+c)}$ 則是表示共有的OTU占總OTU中的比例孽查，其意義則是共有OTU豐度占總OTU豐度的比例，即相似度的一種度量坦喘，而距離=1-相似度盲再。
比如有兩個樣本A和B，OTU及其對應(yīng)豐度信息如表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$

加權(quán)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})$ 本質(zhì)上是對存在于兩個樣本中OTU的進行求和斩例，乘以2的目的是因為兩個樣本的共有的OTU在每個樣本中各自計算1次，分母 $\sum S_{A,i}+\sum S_{B,i}$ 則是兩個樣本中所有OTU個數(shù)之和从橘。

2.1和2.2中的兩個公式本質(zhì)上的意義是一致的樱拴，唯一的不同就是第一個公式是基于OTU的存在與否，而第二個公式是基于OTU的豐度來計算兩個群落微生物的組成差異洋满。

3 UniFrac距離

Bray&Curtis等距離的計算方法往往基于一個前提，就是假設(shè)屬性之間是獨立的珍坊，而事實上往往這些屬性是存在一定關(guān)系的牺勾，因此計算距離的時候也可以把這些屬性之間的關(guān)系考慮在內(nèi)。在宏基因組的研究中阵漏，研究者可以選擇用UniFrac距離來兼顧物種之間的進化關(guān)系計算群落間的距離驻民。
UniFrac距離算法翻具，誕生于2005年底，由Lozupone和Knight開發(fā)[4]回还。
我們可以顧名思義地猜測一下兩位大神作者的取UniFrac這個名字的意圖裆泳。
Uni，即Unique柠硕，獨有的工禾；Frac，即Fraction蝗柔，分?jǐn)?shù)闻葵，也可以看成是比例。
那么具體是如何Uni癣丧，又是怎么個Frac呢槽畔？
請諸位大哥繼續(xù)往下看。
UniFrac距離包括非加權(quán)（Unweighted）和加權(quán)（Weighted）UniFrac距離胁编，這里只介紹非加權(quán)UniFrac距離的原理和計算方法厢钧。

3.1 非加權(quán)UniFrac距離

由于UniFrac距離考慮物種進化關(guān)系，因此在計算UniFrac距離之前嬉橙，我們要先計算其OTU之間的系統(tǒng)發(fā)育關(guān)系早直。
如圖2所示，已知兩個微生物群落樣本憎夷，分別由紅色方塊和綠色圓形表示莽鸿。紫色的分支表示來自某單一群落的后代。UniFrac用紫色樹枝的長度除以樹的總樹枝長度來衡量這兩個群落之間的進化差異拾给。左邊的例子顯示了兩個群落的進化差異很小祥得，因此UniFrac距離相對較低。相反蒋得，右邊的例子顯示了兩個最大差異的群落级及，因為它們最低的共同祖先是樹的根[3]。

圖2 Unweighted UniFrac距離計算方法圖解

非加權(quán)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$ 是在樹上的節(jié)點數(shù)量， $l_i$ 是節(jié)點 $i$ 及其父節(jié)點之間的長度窍侧， $A_i$ 和 $B_i$ 取值 $0$ 或 $1$ 分別表示節(jié)點 $i$ 的子節(jié)點是否存在群落 $A$ 和群落 $B$ 县踢。

UniFrac可以應(yīng)用于多個群落的距離計算，具體方法是：先構(gòu)建一個由所有序列組成的樹伟件，任何一對群落之間的距離則通過修剪掉其他群落的樹枝計算硼啤，這樣樹就只包含來自這對群落的序列。（這在上面給出的公式中是隱式的斧账，因為沒有出現(xiàn)在任何一對群落中的節(jié)點的統(tǒng)計量為零谴返，所以不影響結(jié)果煞肾。）

接下來應(yīng)用隨機排列檢驗計算兩個群落之間的距離是否顯著大于隨機預(yù)期。如圖3所示嗓袱，這是通過將群落標(biāo)簽分配到序列中籍救，同時保持樹不變來實現(xiàn)的。假設(shè)檢驗獲得的 $p$ 值是導(dǎo)致UniFrac距離大于或等于觀測數(shù)據(jù)的隨機排列的概率（圖3）渠抹。

圖3 隨機排列檢驗檢驗群落距離顯著性

兩個群落的距離是這么計算的蝙昙，那么三個的呢？
上述僅介紹了兩個群落之間的距離計算方法逼肯，圖4較為清晰地展示了三個群落計算非加權(quán)UniFrac的具體操作流程耸黑，即A和B構(gòu)建系統(tǒng)發(fā)育樹計算UniFrac距離，C為隨機排列檢驗篮幢，D則是距離矩陣的計算大刊。

圖4 計算Unweighted UniFrac距離矩陣的基本流程。正方形三椿、三角形和圓形表示來自不同群落的序列缺菌。節(jié)點分支如果全部處于單一群落，則為黑色搜锰；如果是節(jié)點分支是多群落共享共享的伴郁，則為灰色。

3.2 加權(quán)UniFrac距離的計算

加權(quán)UniFrac算法實際上可以看作是非加權(quán)UniFrac的擴展版本蛋叼，進一步引入了屬性的量化信息焊傅。
下圖展示了將加權(quán)UniFrac應(yīng)用于兩個由紅色方塊和綠色圓圈表示的的兩個群落樣本的計算結(jié)果。對各個分枝進行加權(quán)狈涮，并計算各群落中分枝下的序列的相對比例狐胎。紫色的分支是那些導(dǎo)致分子非零值的分支「桠桑考慮到序列之間的演化速率差異握巢，加權(quán)UniFrac度量被標(biāo)準(zhǔn)化，使得每個序列對計算的距離貢獻(xiàn)相等松却。

圖5 加權(quán)UniFrac距離的計算

上圖的公式中暴浦， $N$ 是在樹上的節(jié)點數(shù)量， $S$ 是樹所表示的序列的數(shù)目晓锻， $l_i$ 是節(jié)點 $i$ 及其父節(jié)點之間的長度歌焦， $Lj$ 是從樹根到樹梢即序列 $j$ 的總分支長度， $A_i$ 和 $B_i$ 是分別來自群落 $A$ 和 $B$ 的節(jié)點數(shù)量砚哆，而 $A_T$ 和 $B_T$ 分別是群落樣本 $A$ 和 $B$ 序列的總數(shù)独撇。

加權(quán)的UniFrac可以應(yīng)用于多個群落樣本距離矩陣的計算，構(gòu)建一個由所有序列組成的樹。每對群落樣本之間的距離是通過修剪樹來計算的券勺，這樣樹就只包含來自這對群落的序列，可以通過隨機排列測試來測試兩個社區(qū)之間的距離是否大于預(yù)期灿里，這與未加權(quán)的UniFrac類似关炼。

3.2 加權(quán)和非加權(quán)的抉擇

該算法作者在2006年的一篇文章中[6]告訴我們對數(shù)據(jù)集同時進行加權(quán)UniFrac和非加權(quán)UniFrac算法往往會導(dǎo)致不同的結(jié)果，因為它們會導(dǎo)致關(guān)于微生物群落相似性和驅(qū)動微生物多樣性的因素的不同結(jié)論匣吊。
一般來說儒拂，如果群落結(jié)構(gòu)物種的種類差別較大的情況下，用非加權(quán)UniFrac距離要比用非加權(quán)UniFrac距離更能找到區(qū)別色鸳；反之如果僅僅是物種豐度差別大的話社痛，則選擇加權(quán)Unifrac距離。

3.3 Fast UniFrac算法

2010年命雀，Hamady等[7]在Unifrac算法的基礎(chǔ)上開發(fā)了 Fast UniFrac算法蒜哀。如圖4所示，以四個群落樣本（R吏砂、B撵儿、O和G）為例，
1）首先狐血，將各個群落樣本進行根據(jù)每個OTU的有無進行編碼淀歇，比如0號OTU僅僅存在于R中，則編碼為1000匈织；1號OTU僅存在于B中浪默，則編碼為0100，……以此類推缀匕。
2）接下來纳决，去掉暫時不相關(guān)的群落，選取要比較的群落和對應(yīng)的編碼弦追。
3）使用并集算法對環(huán)境集進行比較岳链，將狀態(tài)分配給每個內(nèi)部節(jié)點，比如5號和6號OTU分別編碼為00和01劲件，那么它們共同的上一級父節(jié)點11號節(jié)點則編碼為01掸哑；
4）最后，計算各個群落樣本的枝長百分比再求和零远，就是非加權(quán)Unifrac距離苗分，實際上圖4h中的公式與3.1 非加權(quán)Unifrac距離是等效的。

圖4 UniFrac和Fast UniFrac在流程上的差異（非加權(quán)）牵辣。UniFrac算法：a）環(huán)境以集合的形式存儲在樹對象中摔癣；b）這棵樹被修剪，只包括那些通向想要的環(huán)境的枝干；c）使用集合算法對環(huán)境集進行比較择浊，將狀態(tài)分配給每個內(nèi)部節(jié)點戴卜；d）結(jié)果由另一個樹遍歷計算。Fast UniFrac算法：e）環(huán)境被存儲為提示x環(huán)境計數(shù)的數(shù)組琢岩；f）通過對該數(shù)組進行切片來選擇所選環(huán)境地粪；g）使用數(shù)組片上的數(shù)組操作計算內(nèi)部狀態(tài)魏身；h）將關(guān)聯(lián)數(shù)組和連接兩個環(huán)境中的一個或兩個環(huán)境的節(jié)點的分支長度的乘積相加，從而計算UniFrac值∨榇猓基于數(shù)組的方法可以大大提高效率净赴。

4 寫在后面

最后火诸，想說一說一個文學(xué)作品中永恒的主題——愛情该窗。
愛情在文學(xué)作品中，可以說是突破性非常強的一個東西拌消。
愛情可以突破社會階層挑豌，突破地域，突破宗教信仰拼坎，突破年齡……
作為生物狗浮毯，到這里我都還可以理解。
不過也有我無法理解泰鸡，因為他們跨越了物種债蓝。
很明顯，這不科學(xué)盛龄。

若干年前饰迹，有一首歌火遍大江南北，一首讓我乍一聽就滿臉雞皮疙瘩的歌余舶。我一說名字你準(zhǔn)知道——《狼愛上羊》啊鸭，內(nèi)容通篇如歌名，沒聽過的可以搜一下匿值。

你要說馬和驢之間的愛情赠制，我還可以勉強接受。
雖然染色體數(shù)量不一樣挟憔，無法產(chǎn)生可育后代即騾子钟些，但他們畢竟長相差別不大，也難保馬或驢不臉盲啊绊谭。
狼和羊本是食物鏈關(guān)系政恍，他們是不會有結(jié)果的。
不過达传，還有更過分的篙耗。

曾經(jīng)迫筑，我的某室友QQ昵稱叫“飛鳥和魚”。
啥意思宗弯？這令我很費解脯燃，如果食物鏈，那么為啥不叫“水鳥和魚”蒙保？
后來才知道曲伊，原來是有這么一首詩《飛鳥和魚》，詩里面有這么一句：
“世界上最遙遠(yuǎn)的距離追他，不是我站在你面前，你卻不知道我愛你岛蚤，
而是明明知道彼此相愛邑狸，卻不能在一起〉佣剩”通篇都是這類句型单雾。

又是一對食物鏈的愛情，從生物分類學(xué)角度來講她紫，狼和羊尚同屬哺乳綱硅堆，飛鳥和魚分屬魚綱和鳥綱，已經(jīng)跨綱了……然而贿讹，跨綱算個屁渐逃，還有跨界的。你聽（此處應(yīng)有歌聲）：

“……我愛你愛著你就像老鼠愛大米……”
這歌一定聽過吧民褂，是不是很熟悉茄菊？
又是食物鏈的愛情……。

人類社會中尚且不允許兄妹之間產(chǎn)生愛情赊堪，何況不同物種面殖？
如此，可以總結(jié)一下哭廉，突破了物種的愛情是不會有結(jié)果的脊僚。
我們且看演化生物學(xué)家恩斯特·麥爾的定義：“物種是能夠相互配育的自然種群的類群，這些類群與其它這樣的類群在生殖上相互隔離著遵绰×苫希”

這些距離，由遠(yuǎn)及近街立，總結(jié)起來很有層次感舶衬。
最遠(yuǎn)的，群落層面赎离，Beta多樣性逛犹；
近一點呢，物種間層面，物種和物種之間的距離虽画，不同物種之間性交不可育舞蔽；
再近一點，物種內(nèi)層面码撰，即物種間各個個體之間的距離渗柿，物種內(nèi)性交可育。
還有番外篇脖岛，個體內(nèi)朵栖，基因與基因間的按重組率遺傳距離，以厘摩（cM）為單位柴梆，做過遺傳圖譜的大哥們都知道陨溅。
以上這些均可通過DNA序列相關(guān)數(shù)據(jù)計算獲得。

一般來說绍在，兩個對象之間的距離很直觀门扇，但是意義極為有限，甚至可以說沒有意義偿渡。為了讓這個距離顯得有意義臼寄，現(xiàn)實的分析過程中，往往需要針對多個對象溜宽，分別計算每兩個對象之間的距離吉拳，最后得到一個對稱的距離矩陣（如表4）。盡管應(yīng)用該矩陣可以得知兩個對象之間的距離适揉，但是非常不直觀合武，因而難以對這些對象進行分類。至于如何更直觀地對距離進行可視化展示涡扼，我們就留到后面再說吧稼跳。
表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日發(fā)表于e媛微生態(tài)公眾號。

參考文獻(xiàn)
[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多樣性（二）