微生物群落研究的主要目標是比較不同群落的組成(β多樣性)档痪。在第6章介紹了β多樣性，并舉例說明了如何計算β多樣性指數(shù)邢滑。得到β多樣性指數(shù)后腐螟，可以對其進行統(tǒng)計分析。微生物群研究中的β多樣性分析分為兩類：探索性技術(shù)和有意義的統(tǒng)計檢驗困后。我們在第7章中說明了聚類和排序乐纸。在本章中，我們重點介紹關(guān)于β多樣性的有意義的統(tǒng)計檢驗摇予。已經(jīng)開發(fā)了幾種方法或模型來測試微生物群落組成的差異汽绢。我們將介紹使用置換多變量方差分析(PERMANOVA)、Mantel檢驗侧戴、相似性分析(ANOSIM)宁昭、多響應置換過程(MRPP)和廣義UniFrac距離通過PERMANOVA進行β多樣性的統(tǒng)計檢驗。

9.1 Hypothesis Testing Among Groups Using Permutational Multivariate Analysis of Variance? (PERMANOVA)【使用置換多變量方差分析(PERMANOVA)進行組間假設檢驗】

① PERMANOVA簡介

傳統(tǒng)的多元方差分析(MANOVA)只是一個含有多個因變量的方差分析(ANOVA)酗宋。ANOVA用于檢驗零假設：兩個或多個組之間的均值沒有差異积仗；比較起來，Manova可以用于檢驗零假設：兩個或多個均值向量沒有差異本缠。因此斥扛，它對于多變量數(shù)據(jù)的分析特別強大。然而,傳統(tǒng)的MANOVA假設因變量應該群體內(nèi)部的正態(tài)分布(正態(tài)分布),在所有對因變量的線性關(guān)系,協(xié)變量的全對,在每一個細胞和所有依賴variable-covariate對(線性),因變量展覽同等水平的方差在預測變量的范圍(方差的同質(zhì)性)和方差和協(xié)方差的同質(zhì)性丹锹。由于這些嚴格的假設稀颁，傳統(tǒng)的MANOVA不適合大多數(shù)微生物組多變量數(shù)據(jù)集。例如楣黍，不適合分析微生物組組成與環(huán)境因素(即不同處理組或條件)之間的關(guān)系匾灶。

2001年，Anderson(2001)提出了一種基于平方距離和的分析和劃分的非參數(shù)方法來檢驗兩組或更多組樣本之間沒有差異的假設租漂。它被稱為置換Manova(以前稱為“非參數(shù)Manova”阶女，縮寫為NP-Manova)。Anderson的置換Manova是由相異矩陣表示的哩治。讓我們考慮每對觀測之間的距離矩陣秃踩，設N=an，觀測(點)的總數(shù)业筏，并且dij觀測之間的距離i=1憔杨，…，N蒜胖，并且觀測j=1消别，…抛蚤，N，那么平方和的總和如下所示

總和將距離矩陣(但不包括對角線)的半對角線(或上對角線)中所有距離的平方相加寻狂，然后除以N岁经。SST用于計算所有樣本之間的平均距離。類似地蛇券，組內(nèi)平方和或剩余平方和為

其中ij是指示符缀壤，如果觀測i和觀測j在同一組中，則取值1怀读，否則取值0诉位。也就是說，將同一組中出現(xiàn)的兩個觀測值之間的所有距離的平方相加菜枷，然后除以n苍糠，即每個組中的觀測值數(shù)量。在此基礎上啤誊，利用SSW來計算組內(nèi)樣本間的平均距離岳瞭。用于計算組間平均距離的平方和可以通過用SSW減去SST得到。也就是說蚊锹，SSA=SST-SSW瞳筏。給出了檢驗多元假設的偽F比：

該檢驗統(tǒng)計的基本原理在于，如果來自不同組的點在多變量空間中具有不同的中心位置(在歐幾里德距離的情況下為質(zhì)心)牡昆，與組內(nèi)距離相比姚炕，則組間距離將相對較大，產(chǎn)生的偽F比將相對較大丢烘。從作為F檢驗統(tǒng)計量的公式(9.3)中柱宦，我們可以看到，SSA/(a-1)與SSW/(N-1)之比(稱為信噪比)越大播瞳，F(xiàn)值越大掸刊，從而導致較小的p值。問題是：我們?nèi)绾尾拍塬@得p值赢乓？微生物群和生態(tài)學中的個體變量通常不是正態(tài)分布的忧侧，我們預計歐幾里德距離不一定會用于分析。正如Anderson(2001)注意到的那樣牌芋，即使每個變量都是正態(tài)分布的蚓炬，并且使用了歐幾里德距離，為多變量數(shù)據(jù)計算的平均平方也不會每個都由獨立的v2變量之和組成躺屁，因為盡管單個觀測被認為是獨立的试吁，但是單個物種(OTU或分類群)變量并不是彼此獨立的。因此楼咳，不能使用傳統(tǒng)的表格p值熄捍。Anderson建議對數(shù)據(jù)集中的變量進行隨機打亂(置換)，以生成經(jīng)驗分布母怜。其背后的算法是余耽，在零假設下，各組并沒有真正的不同苹熏，然后多變量觀測(行)將在不同的組之間交換碟贾。可以針對行的所有可能的重新排序重復隨機洗牌轨域，例如1000次袱耽，然后生成特定于該數(shù)據(jù)類型的1000個分布。表示組間顯著性的p值將從經(jīng)驗F分布中獲得「煞ⅲ現(xiàn)在朱巨，讓我們解釋一下如何從排列測試中獲得p值。如上所述枉长，排列意味著將樣本觀測隨機分配到組中冀续。P值是通過將置換的F比與觀察到的F比進行比較來計算的。顯著性檢驗簡單地說就是排列后的F比大于觀察到的F比的分數(shù)必峰。簡單地說洪唐，我們以單向測試為例來解釋排列的實際作用。在單向測試中吼蚁，我們感興趣的是看一個統(tǒng)計數(shù)字是小于還是大于可隨機預期的值凭需。根據(jù)大于或等于觀察到的F統(tǒng)計量的置換偽F統(tǒng)計量的比例來計算p值。換句話說肝匆，我們想知道PERMANOVA之后的置換數(shù)據(jù)集是否產(chǎn)生了相對于實際數(shù)據(jù)集更好的組分辨率粒蜈。如果超過5%的置換F統(tǒng)計量的值大于觀察到的F統(tǒng)計量的值，則p值大于0.05术唬。然后薪伏，我們可以得出結(jié)論，各組之間的任何差異在統(tǒng)計上都不顯著粗仓。與傳統(tǒng)的Manova相比嫁怀，PERMANOVA至少有兩個優(yōu)勢。首先借浊，它不需要任何關(guān)于分布的假設塘淑。其次，它可以使用由任何距離度量計算的距離矩陣蚂斤。顯然翘簇，非分布假設是一種優(yōu)勢吆倦，但后面的特征也很重要伍宦，因為基于歐幾里得距離的分析是計算組內(nèi)樣本之間的平均距離靖秩。換句話說，它是測量群在歐幾里得空間中的中心位置福侈，稱為質(zhì)心。然而，對于許多距離測量兼吓，很難計算出中心位置。例如森枪，在生態(tài)學和微生物學研究中视搏，存在許多情況下半度量的Bray-Curtis測度更合適;然而，在多元Bray-Curtis空間中县袱，我們不能輕易地從樣本數(shù)據(jù)中直接計算出中心位置浑娜。

②用Vegan軟件包實現(xiàn)PERMANOVA

PERMANOVA是通過vegan包中的函數(shù)adonis()實現(xiàn)的。函數(shù)adonis()用于使用半度量和度量距離矩陣分析和劃分平方和式散。adonis()的典型用途包括分析生態(tài)和微生物群落數(shù)據(jù)(按物種(分類群)矩陣采樣)或具有有限數(shù)量的個體樣本和數(shù)千或數(shù)百萬列基因表達的遺傳數(shù)據(jù)筋遭。函數(shù)adonis()允許對由連續(xù)和/或分類預測因子解釋的β多樣性的方差進行類似于方差分析的檢驗。adonis是vegan中推薦的方法杂数。vegan中的其他方法包括MRPP和ANOSIM宛畦。然而，MRPP和ANOSIM都只處理類別預測揍移，并且它們不如adonis好用次和。要實施PERMANOVA，您需要選擇一個距離度量那伐。對于常規(guī)測量踏施，通常選擇“歐幾里得”距離，但對于群落微生物群數(shù)據(jù)罕邀，“Bray”(Bray-Curtis距離)更合適畅形。在文獻中，通常已經(jīng)報告了四個β-分集度量诉探，包括Bray-Curtis距離日熬、Jaccard距離、加權(quán)和未加權(quán)的UniFrac距離肾胯。下面給出一個用法的例子：

adonis(formula, data, permutations = 1000, method = “bray”, contr.unordered =“contr.sum”, contr.ordered = “contr.poly”)?在上述語法中竖席，formula=模型公式，如Y~A+B+C*D敬肚，其中毕荐，Y可以是不同的對象(繼承自“dist”類)、數(shù)據(jù)框或矩陣艳馒；A憎亚、B、C、D可以是因子或連續(xù)變量第美。

在第6章利用VDR小鼠樣本數(shù)據(jù)計算了包括Bray-Curtis指數(shù)在內(nèi)的多種β-多樣性測度蝶锋。在獲得這些β多樣性指數(shù)后，我們可以進行多變量群落分析斋日，以測試微生物群落的組成在不同樣本之間的差異牲览。在這種情況下，我們首先要測試微生物群落的組成是否在糞便和盲腸部位有所不同恶守，并在不同的組(Vdr?/?和WT小鼠)之間存在差異。我們還想在相同的地點(如糞便樣本)測試贡必，與WT組相比兔港，遺傳缺陷組(Vdr?/?小鼠)是否有不同的微生物群組成。

使用Vegan軟件包準備PERMANOVA的實施：因為我們的數(shù)據(jù)集中沒有直接給出位置和組信息仔拟，所以我們需要額外的編程來從樣本ID中提取它們衫樊。以下代碼用于將字符串樣本ID拆分為X1、X2利花、X3三個分量科侈。

> grouping<-data.frame(row.names=rownames(abund_table),t(as.data.frame(strsplit(rownames (abund_table),"_"))))

> grouping<-data.frame(row.names=rownames(abund_table),t(as.data.frame (strsplit(rownames(abund_table),"_"))))

> grouping

X1 X2? ? ? ? X3

5_15_drySt-28F? 5 15 drySt-28F

20_12_CeSt-28F 20 12? CeSt-28F

1_11_drySt-28F? 1 11 drySt-28F

2_12_drySt-28F? 2 12 drySt-28F

3_13_drySt-28F? 3 13 drySt-28F

4_14_drySt-28F? 4 14 drySt-28F

7_22_drySt-28F? 7 22 drySt-28F

8_23_drySt-28F? 8 23 drySt-28F

9_24_drySt-28F? 9 24 drySt-28F

19_11_CeSt-28F 19 11? CeSt-28F

21_13_CeSt-28F 21 13? CeSt-28F

22_14_CeSt-28F 22 14? CeSt-28F

23_15_CeSt-28F 23 15? CeSt-28F

25_22_CeSt-28F 25 22? CeSt-28F

26_23_CeSt-28F 26 23? CeSt-28F

27_24_CeSt-28F 27 24? CeSt-28F

在數(shù)據(jù)集中，“dryST”表示樣品來自糞便炒事，“CeSt”表示樣品來自盲腸臀栈。因此，我們創(chuàng)建了一個位置變量來對來自糞便和盲腸部位的樣本進行分組挠乳。

> grouping$Location <- with(grouping, ifelse(X3%in%"drySt-28F", "Fecal", "Cecal"))

我們進一步從糞便和盲腸部位的WT樣本中分離出Vdr?/?樣本权薯，如下所示：

> grouping$Group <- with(grouping,ifelse(as.factor(X2)%in% c (11,12,13,14,15),c("Vdr-/-"), c("WT")))

檢查變量名并從數(shù)據(jù)集中刪除X1、X2和X3睡扬。

> names(grouping)

> grouping <- grouping[,c(4,5)]

> grouping

Location? Group

5_15_drySt-28F? ? Fecal Vdr-/-

20_12_CeSt-28F? ? Cecal Vdr-/-

1_11_drySt-28F? ? Fecal Vdr-/-

2_12_drySt-28F? ? Fecal Vdr-/-

3_13_drySt-28F? ? Fecal Vdr-/-

4_14_drySt-28F? ? Fecal Vdr-/-

7_22_drySt-28F? ? Fecal? ? WT

8_23_drySt-28F? ? Fecal? ? WT

9_24_drySt-28F? ? Fecal? ? WT

19_11_CeSt-28F? ? Cecal Vdr-/-

21_13_CeSt-28F? ? Cecal Vdr-/-

22_14_CeSt-28F? ? Cecal Vdr-/-

23_15_CeSt-28F? ? Cecal Vdr-/-

25_22_CeSt-28F? ? Cecal? ? WT

26_23_CeSt-28F? ? Cecal? ? WT

27_24_CeSt-28F? ? Cecal? ? WT

在創(chuàng)建位置和組變量之后盟蚣，現(xiàn)在可以調(diào)用函數(shù)adonis()來執(zhí)行PERMANOVA。函數(shù)adonis()的輸入數(shù)據(jù)可以是相異點或數(shù)據(jù)框卖怜；在后一種情況下屎开，adonis()使用vegdist()來查找相異點。為了說明這一點马靠，我們選擇了Bray-Curtis奄抽、Jaccard和S?rensenβ-多樣性度量。他們是在第6章中計算出來的虑粥，可以直接使用如孝。

基因型間Bray-Curtis差異的檢驗差異：實現(xiàn)PERMANOVA的最簡單的R代碼用來檢驗vegan中Vdr?/?和野生型樣本的Bray-Curtis不同之處，如下所示娩贷。

> set.seed(123)

> # adonis="analysis of dissimilarity"?

> adonis(bray ~ Group,data=grouping,permutations = 1000)

Call:

adonis(formula = bray ~ Group, data = grouping, permutations = 1000)

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

Df SumsOfSqs MeanSqs F.Model? ? R2 Pr(>F)

Group? ? ? 1? ? ? 0.23? 0.230? ? 1.5 0.097? 0.18

Residuals 14? ? ? 2.14? 0.153? ? ? ? 0.903? ? ?

Total? ? 15? ? ? 2.37? ? ? ? ? ? ? ? 1.000?

如果沒有給出Bray-Curtis相異度第晰，下面的代碼會給出相同的結(jié)果。

> adonis(abund_table ~ Group,data=grouping,permutations = 1000, method ="bray")

Call:

adonis(formula = abund_table ~ Group, data = grouping, permutations = 1000,? ? ?

method = "bray")

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ?Df? ?SumsOfSqs? ?MeanSqs F.Model? ? R2? ? ?Pr(>F)

Group? ? ? ?1? ? ? 0.23? ? ? ? ? ? ? ? 0.230? ? ? 1.5? ? ? ? ?0.097? ? 0.19

Residuals 14? ? ?2.14? ? ? ? ? ? ? ? 0.153? ? ? ? ? ? ? ? ? ?0.903? ? ?

Total? ? ? ? ? 15? ? ?2.37? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1.000?

除了由于四舍五入導致的p值差異外，這兩個輸出中的所有統(tǒng)計數(shù)據(jù)都是相同的茁瘦。通過排列試驗得到p值品抽。 分組R2為0.097，可用組解釋總方差的9.7%甜熔。對Bray-Curtis距離矩陣的行和列執(zhí)行1000次隨機化圆恤，在零假設下生成R2。在這1,000個值中腔稀，有180個比觀測值0.097大盆昙，因此獲得與觀測值一樣大的值的幾率小于180/1000，表明p值為0.18.焊虏。因此淡喜，我們得出結(jié)論，Vdr?/?和野生型樣本之間的Bray-Curtis距離差異往往是偶然的诵闭。

不同基因型間Jaccard差異的檢驗差異：類似地炼团，我們可以使用PERMANOVA來檢驗VDR?/?和野生型樣本之間的 Jaccard相異度的差異。

> adonis(jaccard ~ Group,data=grouping,permutations = 1000)

Call:

adonis(formula = jaccard ~ Group, data = grouping, permutations = 1000)

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ?Df? ? SumsOfSqs MeanSqs? F.Model? ? R2? ? Pr(>F)

Group? ? ? ?1? ? ? 0.33? ? ? ? ? ? ?0.334? ? ? ? ?1.37? ? ? 0.089? ?0.17

Residuals 14? ? ? 3.40? ? ? ? ? ? ?0.243? ? ? ? ? ? ? ? ? ? ?0.911? ? ?

Total? ? ? ? 15? ? ? 3.74? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1.000?

不同基因型間S?rensen差異的檢驗差異：

> adonis(S?rensen ~ Group,data=grouping,permutations = 1000)

Call:

adonis(formula = S?rensen ~ Group, data = grouping, permutations = 1000)

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ? ?Df? ?SumsOfSqs MeanSqs? ? F.Model? ? R2? ? ? Pr(>F)

Group? ? ? ? ?1? ? ? ?0.119? ? ? ? ? ?0.1191? ? ? ? ?1.3? ? ? ?0.085? 0.21

Residuals? 14? ? ? 1.282? ? ? ? ? ?0.0916? ? ? ? ? ? ? ? ? ? ?0.915? ? ?

Total? ? ? ? ? ?15? ? ?1.401? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1.000

PERMANOVA采用Bray-Curtis疏尿、Jaccard和S?rensenβ多樣性測度的相異度分析(“adonis”)進行瘟芝。根據(jù)BrayCurtis距離(R2=0.097，P=0.18)褥琐、Jaccard距離(R2=0.089锌俱，P=0.17)和S?rensen距離(R2=0.085，P=0.21)踩衩，VDR小鼠和WT小鼠之間的距離無顯著性差異嚼鹉。結(jié)果表明，在樣本量較小的情況下驱富，?/?小鼠和WT小鼠在糞便和盲腸樣本組合的情況下锚赤，這三個差異在0.0 5的統(tǒng)計顯著性水平上沒有顯著差異.

地區(qū)間Bray-Curtis相異度的檢驗差異:以下R代碼用于測試糞便和盲腸樣本之間的Bray-Curtis差異是否不同。

> adonis(bray ~ Location,data=grouping,permutations = 1000)

Call:

adonis(formula = bray ~ Location, data = grouping, permutations = 1000)

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ? ?Df? ?SumsOfSqs MeanSqs? F.Model? ? R2? ? ? Pr(>F)? ?

Location? ? ? 1? ? ? 0.533? ? ? ? ? 0.533? ? ? ? 4.06? ? ? ?0.225? ?0.001 ***

Residuals? ?14? ? ?1.840? ? ? ? ? 0.131? ? ? ? ? ? ? ? ? ? ?0.775? ? ? ? ?

Total? ? ? ? ? ?15? ? ? 2.374? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1.000? ? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R2為0.225褐鸥，說明總方差的22.5%可以由地點來解釋线脚。對Bray-Curtis距離矩陣的行和列執(zhí)行1000次隨機化，在零假設下生成R2叫榕。在這1,000個值中浑侥，有一個比觀測值0.225大，因此獲得與觀測值一樣大的值的幾率小于千分之一晰绎，表明p值為0.001寓落。因此，我們可以得出結(jié)論荞下，糞便和盲腸樣本之間的Bray-Curtis距離差異往往不是偶然的伶选。

地區(qū)間Jaccard差異的檢驗差異：以下R代碼用于測試糞便和盲腸樣本之間的Jaccard差異是否不同：

> adonis(jaccard ~ Location,data=grouping,permutations = 1000)

Call:

adonis(formula = jaccard ~ Location, data = grouping, permutations = 1000)

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ?Df? ? SumsOfSqs? MeanSqs? F.Model? ? R2? ? ? ? Pr(>F)?

Location? ? 1? ? ? ?0.64? ? ? ? ? ? ? ?0.645? ? ? ? 2.92? ? ? 0.173? ? ? 0.006 **

Residuals 14? ? ? 3.09? ? ? ? ? ? ? ?0.221? ? ? ? ? ? ? ? ? ? ?0.827? ? ? ? ?

Total? ? ? ? ?15? ? ? 3.74? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1.000? ? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

不同地點之間S?rensen差異性的檢驗差異：下列R代碼用于測試糞便和盲腸樣本之間的S?rensen差異是否不同：

> adonis(S?rensen ~ Location,data=grouping,permutations = 1000)

Call:

adonis(formula = S?rensen ~ Location, data = grouping, permutations = 1000)

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ? ?Df SumsOfSqs MeanSqs? F.Model? ? R2? ? ? ?Pr(>F)? ?

Location? ? ? 1? ? ?0.352? ? ? ? ? ?0.352? ? ? 4.71? ? ? ?0.252? ? 0.001 ***

Residuals? 14? ? 1.049? ? ? ? ? ?0.075? ? ? ? ? ? ? ? ? ? 0.748? ? ? ? ?

Total? ? ? ? ? ?15? ? 1.401? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1.000? ? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Bray-Curtis距離(R2=0.225史飞，P=0.001)、Jaccard距離(Adonis仰税，R2=0.173构资，P=0.006)和S?rensen距離(Adonis，R2=0.252陨簇，P=0.001)顯示出顯著的分離吐绵。根據(jù)Bray-Curtis距離、Jaccard距離和S?rensen距離的差異河绽，我們可以得出結(jié)論己单，兩個地點之間的差異在統(tǒng)計上是顯著的，并且大約22.5%葵姥、17.3%和25.2%的差異是由地點來解釋的荷鼠。

Bray-Curtis相異度按組對地點排序的貫序檢驗：在群體檢測中，VDR?/?和WT樣本與糞便或盲腸部位無明顯區(qū)別榔幸。在定位檢測中，糞便和盲腸標本與VDR?/?或WT組無明顯區(qū)別矮嫉。實際上削咆，在給定兩個變量組和位置的情況下，我們可以進行排列方差分析的序貫檢驗蠢笋。為了說明這一點拨齐，我們使用了Bray-Curtis相異度來進行貫序檢驗。感興趣的讀者可以使用Jaccard和Jaccard距離測量來嘗試他們自己的測試昨寞。

> adonis(bray ~ Group*Location,data=grouping,permutations = 1000)

Call:

adonis(formula = bray ~ Group * Location, data = grouping, permutations = 1000)

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Df SumsOfSqs MeanSqs F.Model? ? R2? ? ? Pr(>F)? ?

Group? ? ? ? ? ? ? ? ? ? 1? ? 0.230? ? ? ? ? ?0.230? ? ? 1.77? ? ? ? 0.097? ?0.123? ?

Location? ? ? ? ? ? ? ? 1? ? 0.533? ? ? ? ? ?0.533? ? ? ?4.10? ? ? ?0.225? ?0.001 ***

Group:Location? ? ?1? ? 0.051? ? ? ? ? ?0.051? ? ?0.39? ? ? ? 0.021? ? 0.928? ?

Residuals? ? ? ? ? ? ?12? ? 1.559? ? ? ? ?0.130? ? ? ? ? ? ? ? ? ? ? 0.657? ? ? ? ?

Total? ? ? ? ? ? ? ? ? ? ?15? ? 2.374? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1.000? ? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

分組Bray-Curtis相異排序位置的貫序檢驗：

> adonis(bray ~ Location*Group,data=grouping,permutations = 1000)

Call:

adonis(formula = bray ~ Location * Group, data = grouping, permutations = 1000)

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ? ? ? ? ? ?Df? ?SumsOfSqs MeanSqs F.Model? ? R2? ? ? ?Pr(>F)?

Location? ? ? ? ? ? 1? ? ? ? 0.533? ? ? ? ? ?0.533? ? ? ?4.10? ? ? 0.225? ? 0.002 **

Group? ? ? ? ? ? ? ? 1? ? ? ?0.230? ? ? ? ? ?0.230? ? ? ? 1.77? ? ? 0.097? ? 0.095 .

Location:Group? 1? ? ? 0.051? ? ? ? ? ?0.051? ? ? ? 0.39? ? ?0.021? ? ?0.909?

Residuals? ? ? ? ? 12? ? ?1.559? ? ? ? ? ?0.130? ? ? ? ? ? ? ? ? ? ?0.657? ? ? ? ?

Total? ? ? ? ? ? ? ? ? 15? ? ? 2.374? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1.000? ? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

請注意瞻惋，變量的順序很重要。但是援岩，無論是組內(nèi)先試還是先定位歼狼，位置均有統(tǒng)計學意義，p<0.05享怀。

四個水平組間Bray-Curtis差異性的檢驗差異：實際上羽峰，我們可以使用PERMANOVA來檢驗兩種基因型和兩個位點組合的四個水平的群體之間的差異。下面的R代碼創(chuàng)建了一個組合了基因型變量和位置變量的4級組變量添瓷。默認分隔符是‘.’梅屉。

> grouping$Group4<- with(grouping, interaction(Location,Group))

> grouping

? ? ? ? ? ? ? ? ? ? ? ? ? ? Location? Group? ? ? Group4

5_15_drySt-28F? ? Fecal Vdr-/- Fecal.Vdr-/-

20_12_CeSt-28F? ? Cecal Vdr-/- Cecal.Vdr-/-

1_11_drySt-28F? ? Fecal Vdr-/- Fecal.Vdr-/-

2_12_drySt-28F? ? Fecal Vdr-/- Fecal.Vdr-/-

3_13_drySt-28F? ? Fecal Vdr-/- Fecal.Vdr-/-

4_14_drySt-28F? ? Fecal Vdr-/- Fecal.Vdr-/-

7_22_drySt-28F? ? Fecal? ? WT? ? Fecal.WT

8_23_drySt-28F? ? Fecal? ? WT? ? Fecal.WT

9_24_drySt-28F? ? Fecal? ? WT? ? Fecal.WT

19_11_CeSt-28F? ? Cecal Vdr-/- Cecal.Vdr-/-

21_13_CeSt-28F? ? Cecal Vdr-/- Cecal.Vdr-/-

22_14_CeSt-28F? ? Cecal Vdr-/- Cecal.Vdr-/-

23_15_CeSt-28F? ? Cecal Vdr-/- Cecal.Vdr-/-

25_22_CeSt-28F? ? Cecal? ? WT? ? Cecal.WT

26_23_CeSt-28F? ? Cecal? ? WT? ? Cecal.WT

27_24_CeSt-28F? ? Cecal? ? WT? ? Cecal.WT

運用PERMANOVA檢驗四個水平組之間的差異性。

> adonis(bray ~ Group4,data=grouping,permutations = 1000)

Call:

adonis(formula = bray ~ Group4, data = grouping, permutations = 1000)

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ? ? Df SumsOfSqs MeanSqs F.Model? ? R2 Pr(>F)?

Group4? ? ? ? 3? ? 0.814? ? ? ? ? ? ?0.271? ? 2.09? ? ? ? ?0.343? 0.018 *

Residuals? ?12? ? 1.559? ? ? ? ? ? ?0.130? ? ? ? ? ? ? ? ? ? ? 0.657? ? ? ?

Total? ? ? ? ? ?15? ? 2.374? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1.000? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

按Bray-Curtis距離計算鳞贷，樣本分離顯著(adonis坯汤，R2=0.343，P=0.018)搀愧。我們可以得出結(jié)論惰聂，四個水平的組之間的差異在統(tǒng)計學上是顯著的疆偿；根據(jù)Bray-Curtis的不同，大約34.3%的“方差”是由組的差異來解釋的庶近。

Change the Default Dummy Contrasts in R to the Default “Sum” Contrasts in?Vegan Package（ 將R中的默認虛擬對比度更改為Vegan包中的默認“Sum”對比度）：對于設計矩陣中的無序因素翁脆，R默認使用虛擬對比或處理對比。vegan包默認使用“SUM”或ANOVA對比鼻种。無序因子的“總和”對比度用于更改R中的默認對比度反番。

> adonis(bray ~ Group4,data=grouping,permutations = 1000,contr.unordered ="contr.sum")

Call:

adonis(formula = bray ~ Group4, data = grouping, permutations = 1000,

contr.unordered = "contr.sum")

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ? Df SumsOfSqs MeanSqs F.Model? ? R2? ? ? Pr(>F)?

Group4? ? ? 3? ? 0.814? ? ? ? ? 0.271? ? ? ? 2.09? ? ? ? 0.343? 0.009 **

Residuals 12? ? 1.559? ? ? ? ?0.130? ? ? ? ? ? ? ? ? ? ? ?0.657? ? ? ? ?

Total? ? ? ? ?15? ? 2.374? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1.000? ? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

添加有序因素的對比：下面的R代碼為測試添加了有序因子的“poly”對比

> adonis(bray ~ Group4,data=grouping,permutations = 1000,contr.unordered ="contr.sum",contr.ordered = "contr.poly")

Call:

adonis(formula = bray ~ Group4, data = grouping, permutations = 1000,

contr.unordered = "contr.sum", contr.ordered = "contr.poly")

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ?Df SumsOfSqs MeanSqs F.Model? ? R2? ? Pr(>F)?

Group4? ? ? 3? ? 0.814? ? ? ? ? ?0.271? ? ? 2.09? ? ? ?0.343? 0.011 *

Residuals 12? ? 1.559? ? ? ? ? 0.130? ? ? ? ? ? ? ? ? ?0.657? ? ? ?

Total? ? ? ? ?15? ? 2.374? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1.000? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

測試不同群體之間的全局差異：我們可以使用以下兩個示例調(diào)用之一進行置換Manova來測試這四個級別的組之間的全局差異。這兩個調(diào)用都提供相同的輸出叉钥。

> adonis(bray ~ grouping$Group4,permutations = 1000)

Call:

adonis(formula = bray ~ grouping$Group4, permutations = 1000)

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ? ? ? ? ? ? ?Df SumsOfSqs MeanSqs F.Model? ? R2? ? Pr(>F)

grouping$Group4? 3? ? 0.814? ? ? ? ? 0.271? ? ? ? 2.09? ? ? 0.343? 0.017 *

Residuals? ? ? ? ? ? ?12? ? 1.559? ? ? ? ? 0.130? ? ? ? ? ? ? ? ? ? ?0.657? ? ? ?

Total? ? ? ? ? ? ? ? ? ? ?15? ? 2.374? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1.000? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> adonis(abund_table ~ Group4,data=grouping,permutations = 1000, method =

"bray")

Call:

adonis(formula = abund_table ~ Group4, data = grouping, permutations =

1000,method = "bray")

Permutation: free

Number of permutations: 1000

Terms added sequentially (first to last)

? ? ? ? ? ? ? ? ? ? ?Df SumsOfSqs MeanSqs F.Model? ? R2? ?Pr(>F)?

Group4? ? ? ? ? 3? ? 0.814? ? ? ? ? ? 0.271? ? 2.09? ? ? 0.343? 0.01 **

Residuals? ? ?12? ? 1.559? ? ? ? ? 0.130? ? ? ? ? ? ? ? ?0.657? ? ? ? ?

Total? ? ? ? ? ? ?15? ? 2.374? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1.000? ? ? ? ?

---

Signif. codes:? 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

③用RVAideMemoire軟件包實現(xiàn)兩兩置換Manova

在方差分析中罢缸，對置換Manova的顯著檢驗結(jié)果表明，所定義的組之間存在顯著差異投队；但是枫疆，沒有辦法知道哪些組是顯著分開的。在使用vegan進行PERMANOVA之后敷鸦，可以通過使用RVAideMemoire軟件包中的函數(shù)pairwie.perm.manova()來執(zhí)行成對置換Manova息楔，以進行每個組水平的成對比較和多個測試的校正。此函數(shù)的一個調(diào)用示例如下：

pairwise.perm.manova(resp, fact, test = c(“Pillai”, “Wilks”,”Hotelling-Lawley”,“Roy”, “Spherical”), nperm = 1000, progress = TRUE, p.method = ”fdr”)

其中扒披，resp = response值依、典型矩陣(每個變量一列)或距離矩陣或數(shù)據(jù)框；fact = grouping factor碟案；test = choice of test statistic when response is a matrix;愿险；nperm=獲得P值的排列數(shù)目；progress=指示是否應該顯示進度條的邏輯价说；以及p.method=用于校正P值的方法辆亏。

在多重比較中使用方差分析、Kruskal-Wallis Test鳖目。在兩兩排列的Manova中提供了幾種p值調(diào)整方法扮叨，包括“bonferroni” , “holm” ， “hochberg” and “hommel” ,“BH” or? “fdr” , “BY”疑苔。此包中不提供Tukey方法甫匹。Type? P.adjust()檢查R中的選項和引用。如果您不想調(diào)整p值惦费，請使用傳遞選項(“None”)兵迅。要執(zhí)行兩兩排列的Manova，我們首先需要安裝和加載RVAideMemoire包薪贫。

> install.packages("RVAideMemoire")

> library(RVAideMemoire)

然后恍箭，我們可以使用以下三個示例調(diào)用中的任何一個來執(zhí)行兩兩排列的Manova。

> set.seed(0)

> pairwise.perm.manova(bray,grouping$Group4,nperm=1000)

> # o r

> pairwise.perm.manova(vegdist(abund_table,"bray"),grouping$Group4,nperm=1000)

> # o r

> pairwise.perm.manova(bray, grouping$Group4, test = c("Pillai","Wilks","Hotelling-Lawley", "Roy", "Spherical"),nperm = 1000,

+? progress = TRUE, p.method = "fdr")

使用“none“方法校正P值：當使用p.method=“None”時瞧省，p值不會校正扯夭。

> pairwise.perm.manova(bray, grouping$Group4, test = c("Pillai", "Wilks",

"Hotelling-Lawley", "Roy", "Spherical"), nperm = 1000,

+? ? ? ? ? ? ? ? ? ? ? progress = TRUE, p.method = "none")

Pairwise comparisons using permutation MANOVAs on a distance matrix

data:? bray by grouping$Group4

1000 permutations

? ? ? ? ? ? ? ? ? ? Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.02? ? ? ? ? ? ? ? ? ?-? ? ? ? ? ? ? ? ? ? -

Cecal.WT? ? 0.64? ? ? ? ? ? ? ? ? ? 0.07? ? ? ? ? ? ? -

Fecal.WT? ? 0.08? ? ? ? ? ? ? ? ? ?0.16? ? ? ? 0.50? ?

P value adjustment method: none

用“Bonferroni”法校正P值

> pairwise.perm.manova(bray, grouping$Group4, test = c("Pillai", "Wilks","Hotelling-Lawley", "Roy", "Spherical"), nperm = 1000,

+? ? ? ? ? ? ? ? ? ? ? progress = TRUE, p.method = "bonferroni")

Pairwise comparisons using permutation MANOVAs on a distance matrix

data:? bray by grouping$Group4

1000 permutations

? ? ? ? ? ? ? ? ? ? ?Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.1? ? ? ? ? ? ? ? ? ? ? ? -? ? ? ? ? ? ? ? ? -

Cecal.WT? ? 1.0? ? ? ? ? ? ? ? ? ? ? 0.3? ? ? ? ? ? ? ? -

Fecal.WT? ? 0.4? ? ? ? ? ? ? ? ? ? 1.0? ? ? ? ? 1.0? ?

P value adjustment method: bonferroni

用“holm”法校正P值

> pairwise.perm.manova(bray, grouping$Group4, test = c("Pillai", "Wilks","Hotelling-Lawley", "Roy", "Spherical"), nperm = 1000,

+? ? ? ? ? ? ? ? ? ? ? progress = TRUE, p.method = "holm")

Pairwise comparisons using permutation MANOVAs on a distance matrix

data:? bray by grouping$Group4

1000 permutations

? ? ? ? ? ? ? ? ? ? Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.09? ? ? ? ? ? ? ? ? ? ?-? ? ? ? ? ? ? -

Cecal.WT? ? 1.00? ? ? ? ? ? ? ? ? ?0.24? ? ? ? ? -

Fecal.WT? ? 0.27? ? ? ? ? ? ? ? ? 0.52? ? ? ? ?1.00? ?

P value adjustment method: holm

用“Hochberg”法校正P值

> pairwise.perm.manova(bray, grouping$Group4, test = c("Pillai", "Wilks","Hotelling-Lawley", "Roy", "Spherical"), nperm = 1000,

+? ? ? ? ? ? ? ? ? ? ? progress = TRUE, p.method = "hochberg")

Pairwise comparisons using permutation MANOVAs on a distance matrix

data:? bray by grouping$Group4

1000 permutations

? ? ? ? ? ? ? ? ? ?Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.06? ? ? ? ? ? ? ? -? ? ? ? ? ? ? ? ? ? ?-

Cecal.WT? ? 0.65? ? ? ? ? ? ? 0.23? ? ? ? ? ? ? ? -

Fecal.WT? 0.26? ? ? ? ? ? ? ? ? 0.45? ? ? ? ? ?0.65? ?

P value adjustment method: hochberg

用“Hommel”法校正P值

> pairwise.perm.manova(bray, grouping$Group4, test = c("Pillai", "Wilks","Hotelling-Lawley", "Roy", "Spherical"), nperm = 1000,

+? ? ? progress = TRUE, p.method = "hommel")

Pairwise comparisons using permutation MANOVAs on a distance matrix

data:? bray by grouping$Group4

1000 permutations

? ? ? ? ? ? ? ? ? ? ?Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.07? ? ? ? ? ? ? ? ? ?-? ? ? ? ? ? ? ? ? -

Cecal.WT? ? 0.64? ? ? ? ? ? ? ? ?0.20? ? ? ? ? ? ? ? -

Fecal.WT? ? 0.30? ? ? ? ? ? ? ? ?0.53? ? ? ? ? ? 0.64? ?

P value adjustment method: hommel

使用“BH”方法校正P值

> pairwise.perm.manova(bray, grouping$Group4, test = c("Pillai", "Wilks","Hotelling-Lawley", "Roy", "Spherical"), nperm = 1000,

+? ? ? ? ? ? ? ? ? ? ? progress = TRUE, p.method = "BH")

Pairwise comparisons using permutation MANOVAs on a distance matrix

data:? bray by grouping$Group4

1000 permutations

? ? ? ? ? ? ? ? ? ? ?Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.1? ? ? ? ? ? ? ? ? ? ?-? ? ? ? ? ? ? ? ? ? ? ?-

Cecal.WT? ? 0.7? ? ? ? ? ? ? ? ? ? 0.2? ? ? ? ? ? ? ? ? ? -

Fecal.WT? ? 0.2? ? ? ? ? ? ? ? ? ?0.2? ? ? ? ? ? ? ? ? 0.6? ?

P value adjustment method: BH

使用“FDR”方法校正P值

> pairwise.perm.manova(bray, grouping$Group4, test = c("Pillai", "Wilks","Hotelling-Lawley", "Roy", "Spherical"), nperm = 1000,

+? ? ? ? ? ? ? ? ? ? ? progress = TRUE, p.method = "fdr")

Pairwise comparisons using permutation MANOVAs on a distance matrix

data:? bray by grouping$Group4

1000 permutations

? ? ? ? ? ? ? ? ? Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.1? ? ? ? ? ? ? ? ? ? ? ?-? ? ? ? ? ? ? -

Cecal.WT? ? 0.6? ? ? ? ? ? ? ? ? ? ?0.1? ? ? ? ? -

Fecal.WT? ? 0.1? ? ? ? ? ? ? ? ? ? ? ?0.3? ? ? 0.6? ?

P value adjustment method: fdr

使用“BY”方法校正P值

> pairwise.perm.manova(bray, grouping$Group4, test = c("Pillai", "Wilks","Hotelling-Lawley", "Roy", "Spherical"), nperm = 1000,

+? ? ? ? ? ? ? ? ? ? ? progress = TRUE, p.method = "BY")

Pairwise comparisons using permutation MANOVAs on a distance matrix

data:? bray by grouping$Group4

1000 permutations

? ? ? ? ? ? ? ? ? Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.3? ? ? ? ? ? ? ? ? -? ? ? ? ? ? ? ? -

Cecal.WT? ? 1.0? ? ? ? ? ? ? ? ? ? 0.3? ? ? ? ? -

Fecal.WT? ? 0.3? ? ? ? ? ? ? ? ? 0.6? ? ? ? ? 1.0? ?

P value adjustment method: BY

? ? ? ? ? ? ? ? ?Cecal.Vdr-/- Fecal.Vdr-/- Cecal.WT

Fecal.Vdr-/- 0.02? ? ? ? ? ? ? ? ? ? -? ? ? ? ? ? ? ?-

Cecal.WT? ? 0.64? ? ? ? ? ? ? ? ?0.07? ? ? ? ? ? -

Fecal.WT? ? 0.08? ? ? ? ? ? ? ? ?0.16? ? ? ? 0.50? ?

P value adjustment? method: none

未經(jīng)校正鳍贾，糞便VdR?/?和盲腸VdR?/?分離顯著，p值為0.0 2交洗。但是骑科，使用這兩種方法進行調(diào)整后，沒有兩個配對項目在0.05的顯著性水平上顯著构拳。

④Test Group Homogeneities Using the Function betadisper()（使用函數(shù)betadisper()測試組同質(zhì)性）

在使用函adonis()測試組平均值差異之后咆爽，我們可以使用函數(shù)betadisper()來測試組同質(zhì)性的差異。adonis()是類似于多元方差分析置森，betadisper()類似于Levene的方差相等檢驗斗埂。

> homo <-with(grouping,betadisper(bray, Group4))

> homo

Homogeneity of multivariate dispersions

Call: betadisper(d = bray, group = Group4)

No. of Positive Eigenvalues: 11

No. of Negative Eigenvalues: 4

Average distance to median:

Cecal.Vdr-/- Fecal.Vdr-/-? Cecal.WT? ? Fecal.WT

? ? 0.342? ? ? ? 0.285? ? ? ? 0.309? ? ? ? 0.239

Eigenvalues for PCoA axes:

PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8

0.775? ? 0.628? ? 0.358? ? 0.174? ? 0.161? ?0.116? 0.100? ? 0.045

該函數(shù)具有用于圖形顯示的plot和boxplot方法。下面的R代碼產(chǎn)生圖9.1中的PCoA圖凫海。

> plot(homo)

下面的R代碼生成圖9.2中四組到質(zhì)心的距離的箱線圖呛凶。

> boxplot(homo)

我們既可以使用標準參數(shù)方差分析，也可以使用置換檢驗(置換檢驗)來分析擬合模型的顯著性行贪。

> anova(homo)

Analysis of Variance Table

Response: Distances

? ? ? ? ? ? ? ? ? ?Df? Sum Sq Mean Sq? ? ?F value Pr(>F)

Groups? ? ? ?3? ? ?0.021? ? ? 0.00701? ? 0.54? ? ? 0.67

Residuals? 12? ? 0.157? ? ?0.01305? ? ? ? ? ? ?

> permutest(homo)

Permutation test for homogeneity of multivariate dispersions

Permutation: free

Number of permutations: 999

Response: Distances

? ? ? ? ? ? ? ? ? ?Df? ?Sum Sq? ?Mean Sq? ?F? ? N.Perm? Pr(>F)

Groups? ? ? ?3? ? ? 0.021? ? ?0.00701? ?0.54? ? 999? ? 0.67

Residuals? 12? ? ?0.157? ? ?0.01305?

參數(shù)檢驗和排列檢驗均表明漾稀，多元離散度在0.05顯著性水平上無統(tǒng)計學差異。此外建瘫，使用參數(shù)Tukey‘s HSD檢驗可以分析組之間的成對差異县好，如下所示：

> TukeyHSD(homo)

Tukey multiple comparisons of means

95% family-wise confidence level

Fit: aov(formula = distances ~ group, data = df)

$group

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?diff? ? ? lwr? ? ? ? ?upr? ? ?p adj

Fecal.Vdr-/--Cecal.Vdr-/- -0.05622 -0.2707 0.1583 0.8629

Cecal.WT-Cecal.Vdr-/-? -0.03273 -0.2804 0.2150 0.9786

Fecal.WT-Cecal.Vdr-/-? -0.10244 -0.3501 0.1453 0.6221

Cecal.WT-Fecal.Vdr-/-? 0.02349 -0.2242 0.2712 0.9918

Fecal.WT-Fecal.Vdr-/-? -0.04623 -0.2939 0.2015 0.9437

Fecal.WT-Cecal.WT? ? ? ? -0.06972 -0.3466 0.2072 0.8760

結(jié)果表明，四組間的配對比較無統(tǒng)計學差異暖混。因此，我們可以得出結(jié)論翁授，這四組的方差是均勻的拣播。

9.2 Hypothesis Tests Among Group-Differences Using Mantel Test (MANTEL)（使用Mantel檢驗(Mantel)進行組間差異的假設檢驗）

①不同矩陣的Mantel檢驗和部分Mantel檢驗介紹

Mantel發(fā)起了一種置換測試程序來測試兩個距離矩陣之間的相關(guān)性，并在Mantel和Valand中得到進一步發(fā)展收擦。由于這一開端贮配，這一過程在生物和環(huán)境科學中被稱為Mantel檢驗，在統(tǒng)計文獻中也被稱為Mantel和Valand的非參數(shù)Manova塞赂。通常泪勒，相關(guān)分析用于量化兩個連續(xù)變量之間的關(guān)聯(lián)，或者是自變量和因變量之間的關(guān)聯(lián)宴猾，或者是兩個自變量之間的關(guān)聯(lián)圆存。Mantel的檢驗具有優(yōu)勢：它可以應用于不同類型的變量，包括分類仇哆、排名或區(qū)間尺度數(shù)據(jù)沦辙，因為它使用距離(相異)矩陣作為輸入數(shù)據(jù)。Mantel檢驗的設置是回歸分析讹剔，其中變量本身是總結(jié)樣本位置之間成對相似性的距離或相異矩陣油讯。在這種設置下详民，因變量不是“樣本i上分類單元X的豐度”，而可能是“樣本i和j上分類單元X的平均數(shù)量的相似性”陌兑；同樣沈跨，預測變量也不是單個樣本的“條件”，而可能是樣本之間的“條件相似性”兔综《隽荩總而言之，Mante檢驗是一種相關(guān)性分析邻奠，因為它分析了兩個變量或距離矩陣之間的關(guān)聯(lián)笤喳。它僅僅是距離矩陣的回歸。實際上碌宴，Mantel檢驗的威力和通用性在于它處理了距離矩陣杀狡，并構(gòu)建了回歸分析框架。在相關(guān)分析中包含分類變量贰镣，并將這些變量轉(zhuǎn)換為距離(差異)指標的能力呜象，使指標更好地用于假設檢驗，這對生態(tài)學家和微生物研究人員特別有用碑隆。因此恭陡，Mantel的測試被認為克服了解釋物種-環(huán)境關(guān)系或一般分類-環(huán)境關(guān)系的一些固有問題。通常上煤，Mantel檢驗中的零假設是這樣構(gòu)造的：響應變量矩陣中樣本之間的距離與另一個解釋變量矩陣中樣本之間的距離不是線性相關(guān)的休玩。我們可以問的操作性問題是這樣的：“在預測變量(環(huán)境)變量方面相似的樣本在因變量(分類)變量方面也傾向于相似嗎？”或者“相近的樣品在成分上也是相似的嗎劫狠？”或者換一種說法拴疤，“不同類群或環(huán)境不同的樣本，在分類群組成或成分上是否也不同独泞？”測試統(tǒng)計如下：

其中呐矾，z=檢驗統(tǒng)計量，是Hadamard乘積(矩陣)懦砂。在乘積中蜒犯，具有相同維數(shù)的兩個矩陣產(chǎn)生另一個矩陣，其中每個元素ij是原始兩個矩陣的元素ij的乘積荞膘；yi是相依距離矩陣或凝聚向量罚随；xi是分組的預測矩陣或?qū)Ρ染仃嚕绻麡颖驹跇颖窘M中衫画，則為預測矩陣或?qū)Ρ染仃嚭谅蝗绻麡颖驹诓煌慕M中，則為1削罩。適當?shù)木嚯x度量可以是單變量(例如瞄勾，“類桿菌豐度相似性”)费奸，也可以是多變量(例如，Bray-Curtis‘s进陡、Jaccard’s愿阐、S?renson‘s相似性指數(shù))；n是樣本數(shù)量趾疚∮Ю可以將Mantel檢驗過程分為以下四個主要步驟：第一步：計算相異矩陣。第二步：計算檢驗統(tǒng)計量Z糙麦。第三步：檢驗顯著性辛孵。

Mantel檢驗通過判斷一組變量中的相似(或接近)是否與另一組變量中的相似(或接近)相關(guān)颠黎，來比較兩組不同之處柑晒。基本上喜滨，它是不同條目之間的相關(guān)性焚廊。由于距離矩陣的元素不是獨立的冶匹，且N個樣本(N(N?1)/2)之間存在太多差異，因此不適用常規(guī)顯著性檢驗咆瘟。Mantel開發(fā)了漸近測試統(tǒng)計嚼隘，但在素食包中，函數(shù)mantel()使用排列測試袒餐。與PERMANOVA檢驗一樣飞蛹，通過隨機重新排列輸入距離矩陣之一(第一相異度矩陣)的行和列來獲得檢驗統(tǒng)計量z的p值。顯著性檢驗就是排列后的Z大于觀察到的Z的分數(shù)灸眼。它是給定這個大小或更大的Z的概率桩皿。如果隨機化經(jīng)常產(chǎn)生大于或等于觀測數(shù)據(jù)的相關(guān)性，幾乎沒有證據(jù)表明相關(guān)性不同于零幢炸。如果對這兩個距離進行秩變換，則Mantel檢驗與ANOSIM相同拒贱，與秩變換MRPP相似宛徊。

第四步：計算兩個矩陣的相關(guān)系數(shù)r，確定它們之間的關(guān)系強度逻澳。標準化的Mantel統(tǒng)計數(shù)據(jù)r如下所示：

將標準化Mantel統(tǒng)計量r解釋為與其他類型的相關(guān)系數(shù)(如皮爾遜系數(shù)r)一樣是“效應大小”的量度的相關(guān)系數(shù)闸天。r的值落在?1?到+1的范圍內(nèi)，接近?1表示強負相關(guān)斜做，+1表示強正相關(guān)苞氮。如果樣本之間的平均組內(nèi)距離小于樣本之間的平均總距離，則r>0瓤逼；r越大笼吟，兩組之間的相關(guān)性越強库物。Mantel檢驗用于比較兩個距離(不相似)矩陣(例如A和B)之間的相關(guān)性，而部分Mantel檢驗用于估計這兩個矩陣之間的相關(guān)性贷帮，同時控制控制矩陣C的影響戚揭。目標是去除虛假相關(guān)性。在生態(tài)學文獻中撵枢，部分Mantel檢驗的一個典型例子用于將群落距離矩陣與從環(huán)境參數(shù)導出的另一個距離矩陣進行比較民晒，同時使用地理距離作為第三個“控制”距離矩陣。在微生物組研究中锄禽，第三基質(zhì)C可以是環(huán)境基質(zhì)潜必，也可以是由環(huán)境基質(zhì)或另一變量(如處理成員或條件)創(chuàng)建的實驗設計矩陣。但請記住沃但，從環(huán)境參數(shù)導出的環(huán)境矩陣和設計矩陣本質(zhì)上都需要是數(shù)值/連續(xù)的磁滚；還需要是距離(不同)矩陣。如果創(chuàng)建的矩陣不是距離(相異)矩陣绽慈，則在將其用作部分Mantel測試中的控制矩陣之前恨旱，應使用距離函數(shù)將其轉(zhuǎn)換為距離(相異)矩陣。部分Mantel檢驗是通過兩個步驟來構(gòu)造的：首先坝疼，進行A和C之間的回歸以構(gòu)造殘差矩陣A‘搜贤，并進行B和C之間的回歸以構(gòu)造殘差矩陣B’；然后通過標準的Mantel檢驗來比較兩個殘差矩陣A‘和B’钝凶。

②用Vegan軟件包演示Mantel檢驗

兩個群落距離矩陣的相關(guān)性檢驗：隨著深度采樣的進行仪芒，Bray-Curtis、Jaccard和S?rensen相異度指數(shù)通常具有較高的相關(guān)性耕陷。我們可以對這些矩陣進行 mantel檢驗掂名。mantel檢驗可以通過vegan、ape和ade4包進行哟沫。在這里饺蔑，我們使用vegan包進行測試。一種語法如下所示：mantel(xdis, ydis, method = “pearson”, permutations = 1000)

其中嗜诀，xdis猾警、ydis是相異矩陣或距離對象；methods是相關(guān)方法隆敢，是字符串发皿，被接受為“Pearson”、“Spearman”或“Kendall”拂蝎；permutations是在評估顯著性時需要指定的排列數(shù)量穴墅。在第7章。我們使用GUniFrac軟件包中的吸煙者數(shù)據(jù)集繪制了樹，并說明了排序技術(shù)玄货。在本節(jié)中皇钞，我們再次使用它來說明Mantel檢驗。以下R代碼加載包誉结、訪問和子集數(shù)據(jù)鹅士。

> library(GUniFrac)

> data(throat.otu.tab)

> otu_table <-throat.otu.tab

> data(throat.meta)

> data(throat.tree)

> library(dplyr)

> throat_meta <- select(throat.meta, SmokingStatus, Age, Sex, PackYears)

Bray-Curtis和Jaccard相異性的Pearson積矩相關(guān)檢驗：Bray-Curtis和Jaccard相異指數(shù)的相關(guān)性估計如下。

> library(vegan)

> mantel(bray, jaccard,"pearson",permutations=1000)

Mantel statistic based on Pearson's product-moment correlation

Call:

mantel(xdis = bray, ydis = jaccard, method = "pearson", permutations = 1000)

Mantel statistic r: 0.993

Significance: 0.001

Upper quantiles of permutations (null model):

90%? 95% 97.5%? 99%

0.154 0.201 0.233 0.266

Permutation: free

Number of permutations: 1000

Mantel檢驗可以這樣解釋：我們指定排列=1000來對Bray-Curtis距離矩陣的行和列執(zhí)行1000次隨機化惩坑，其在零假設下生成相關(guān)性分布：Bray-Curtis矩陣中的樣本之間的距離與Jaccard矩陣不是線性相關(guān)的掉盅。結(jié)果表明，在這種情況下以舒，在這1,000個值中趾痘，沒有一個大于觀測值0.993(即排列數(shù)<觀察值=1，>排列數(shù)>觀察值=0蔓钟，排列數(shù)=觀察值=1)永票，因此獲得與觀察值一樣大的值的幾率小于千分之一，表示p值為0.001滥沫。由于偶然產(chǎn)生的大關(guān)聯(lián)是如此之小侣集，我們可以得出結(jié)論，Bray-Curtis距離差異隨Jaccard距離線性增加兰绣，換言之世分，這兩個差異指數(shù)呈高度正相關(guān)。我們可以繪制Bray-Curtis和Jaccard矩陣的pearson相關(guān)圖缀辩，如下所示：

> plot(bray, jaccard, main="Scatter plot of Bray-Curtis index vs. Jaccard index")

pearson系數(shù)r為0.993臭埋，p值為0.001。Bray-Curtis和Jaccard矩陣中元素的散點圖表明這兩個距離/不同之間存在線性關(guān)系臀玄，如圖9.3所示瓢阴。

Bray-Curtis和S?rensen異同的Spearman秩相關(guān)檢驗：我們使用Spearman方法計算了Bray-Curtis和S?rensen相異指數(shù)的相關(guān)性，如下所示：

> mantel(bray, S?rensen,"spearman",permutations=1000)

Mantel statistic based on Spearman's rank correlation rho

Call:

mantel(xdis = bray, ydis = S?rensen, method = "spearman", permutations

= 1000)

Mantel statistic r: 0.514

Significance: 0.001

Upper quantiles of permutations (null model):

90% 95% 97.5% 99%

0.132 0.178 0.206 0.241

Permutation: free

Number of permutations: 1000

Jaccard?和?S?rensen異同的Kendall秩相關(guān)檢驗：下面健无，我們使用Kendall方法計算Jaccard和S?rensen相異指數(shù)之間的相關(guān)性荣恐。

> mantel(jaccard,S?rensen,"kendall",permutations=1000)

Mantel statistic based on Kendall's rank correlation tau

Call:

mantel(xdis = jaccard, ydis = S?rensen, method = "kendall", permutations

= 1000)

Mantel statistic r: 0.358

Significance: 0.001

Upper quantiles of permutations (null model):

90%? 95%? 97.5%? ? 99%

0.0961 0.1285 0.1521 0.1849

Permutation: free

Number of permutations: 1000

測試群落矩陣和實驗設計的相關(guān)性：之前Mantel檢驗在兩個群落矩陣之間進行。在這里累贤，我們選擇Bray-Curtis矩陣作為相依距離矩陣募胃，元數(shù)據(jù)中的吸煙狀況作為實驗設計矩陣來說明檢驗。我們將測試吸煙者和不吸煙者的樣本之間的Bray-Curtis不同之處是否不同畦浓。變量SmokingStatus在喉嚨中元數(shù)據(jù)集。以下R代碼提取此變量并將其轉(zhuǎn)換為距離矩陣检疫。

> library(dplyr)

> group <-select(throat.meta, SmokingStatus)

> group$Status <- with(group, ifelse(SmokingStatus%in%"Smoker", 1, 0))

> group <- group[,-1]

> group_dist <- vegdist(scale(group), "euclid")

下面的R代碼對Bray-Curtis矩陣(依賴距離矩陣)和吸煙狀況(預測矩陣)的pearson積矩相關(guān)性進行Mantel檢驗讶请。

> mantel(group_dist,bray,"pearson",permutations = 1000)

Mantel statistic based on Pearson's product-moment correlation

Call:

mantel(xdis = group_dist, ydis = bray, method = "pearson", permutations= 1000)

Mantel statistic r: 0.0728

Significance: 0.003

Upper quantiles of permutations (null model):

90%? ? 95%? 97.5%? ? 99%

0.0245 0.0336 0.0436 0.0581

Permutation: free

Number of permutations: 1000

Mantel統(tǒng)計量r為0.0728，p值為0.003，表明結(jié)果在a為0.05時具有統(tǒng)計學意義夺溢。因此论巍，我們得出結(jié)論，吸煙狀況確實可以預測Bray-Curtis差異风响，盡管這兩個矩陣之間的正相關(guān)性相對較弱嘉汰。通過指定1000個排列來獲得p值。

控制第三矩陣的兩個距離矩陣的相關(guān)性的部分Mantel檢驗：Mantel部分測試用于在控制環(huán)境距離矩陣的同時進行社區(qū)距離矩陣和設計距離矩陣的相關(guān)性分析状勤。例如鞋怀，Mantel部分測試可以用來確定當控制Bray-Curtis距離(相異度)矩陣和吸煙狀況距離矩陣時，該矩陣與年齡持搜、性別和每年的包裝是否有顯著的相關(guān)性密似。假設我們想知道“在控制環(huán)境矩陣(預測矩陣：每年的年齡、性別和包裝)的同時葫盼，Bray-Curtis相異度(依賴矩陣)中有多少變異性是由吸煙狀況(設計矩陣)解釋的残腌？”然后，通過構(gòu)造吸煙狀況與每年年齡贫导、性別和包裝的預測矩陣之間的回歸的殘差矩陣A‘和Bray-Curtis相異度與每年的年齡抛猫、性別和包裝的預測矩陣之間的回歸的殘差矩陣B’，來計算檢驗統(tǒng)計量孩灯，該檢驗統(tǒng)計量是通過構(gòu)造吸煙狀況與每年的年齡闺金、性別和包裝的預測矩陣之間的回歸的殘差矩陣A‘和Bray-Curtis相異度與每年的預測矩陣之間的回歸的殘差矩陣B’來計算的。然后钱反，通過標準Mantel檢驗比較兩個殘差矩陣A‘和B’掖看。下面的R代碼子集是喉部元數(shù)據(jù)集中的環(huán)境數(shù)據(jù)集。

> meta <- select(throat.meta, Age,Sex,PackYears

在原始數(shù)據(jù)集中面哥，性別是一個標記為“Male”和“Female”的字符變量哎壳。為了將此數(shù)據(jù)集轉(zhuǎn)換為數(shù)字距離矩陣，我們需要把它重新編碼成一個數(shù)值尚卫。以下R代碼創(chuàng)建數(shù)字變量性別归榕。

> meta$Gender <- with(meta,ifelse(Sex%in%"Male", 1, 0))

> env <-select(meta, Age,Gender,PackYears)

下面的R代碼使用vegan包中的函數(shù)vegdist()將此數(shù)據(jù)集轉(zhuǎn)換為距離矩陣。

> env_dist <- vegdist(scale(env), "euclid")

Mantel部分測試使用pearson方法進行吱涉，如下所示刹泄。

> mantel.partial(group_dist,bray,env_dist, method = "pearson", permutations=1000)

Partial Mantel statistic based on Pearson's product

Call:

mantel.partial(xdis = group_dist, ydis = bray, zdis = env_dist,? ? method="pearson", permutations = 1000)

Mantel statistic r: 0.068

Significance: 0.008

Upper quantiles of permutations (null model):

90%? ? 95%? 97.5%? ? 99%

0.0215 0.0321 0.0438 0.0641

Permutation: free

Number of permutations: 1000

0.068的Mantel統(tǒng)計表明，在控制環(huán)境矩陣隨年齡怎爵、性別和包次的差異的情況下特石，吸煙狀況與Bray-Curtis距離(不相似)矩陣之間的正相關(guān)關(guān)系相對不強。然而鳖链，p值為0.008姆蘸，表明結(jié)果在a=0.05時具有統(tǒng)計學意義。

9.3 Hypothesis Tests Among-Group Differences Using ANOSIM（假設檢驗組間差異的ANOSIM方法）

①相似度分析(ANOSIM)簡介

ANOSIM只是基于兩個距離矩陣之間的標準化秩相關(guān)的Mantel檢驗的修改版本。ANOSIM檢驗是由發(fā)展起來的逞敷，它是社區(qū)生態(tài)學家和微生物組研究人員經(jīng)常使用的一種無分布的多變量數(shù)據(jù)分析方法狂秦。它是一種基于組內(nèi)和組內(nèi)相似性的排列檢驗來檢驗兩組或更多組樣本之間沒有差異的假設的非參數(shù)程序。它根據(jù)分組因子或?qū)嶒炋幚硭酵凭瑁容^不同采樣單位之間物種(或任何其他分類群)豐度和組成的變化(例如裂问，β多樣性)。與方差分析類似牛柒，ANOSIM將組成員或治療水平視為因素堪簿，并將其建模為解釋變量。相似性分析基于一個簡單的概念：如果被測試組是有意義的焰络，那么組內(nèi)的樣本在組成上應該比來自不同組的樣本更相似戴甩。因此，零假設是：試驗組成員或治療條件之間沒有差異闪彼。本試驗采用Bray-Curtis相似性測度方法甜孤。ANOSIM檢驗統(tǒng)計量基于組間和組內(nèi)平均排名的差異。具體如下：

執(zhí)行ANOSIM主要有五個步驟：第一步：計算相異矩陣畏腕。步驟2：計算等級相異度缴川，并將相異度最小的等級分配為1。步驟3：計算組內(nèi)和組內(nèi)等級差異的平均值描馅。

R被解釋為與其他類型的相關(guān)系數(shù)(如pearson系數(shù))一樣把夸，是“效應大小”的量度的相關(guān)系數(shù)。檢驗統(tǒng)計量是檢驗在零假設下不同組間是否存在差異铭污。如果零假設是正確的恋日，則R=0，這表明組內(nèi)和組內(nèi)的差異在平均上是相同的嘹狞。當高相似度和低相似度完美地混合在一起岂膳，并且與群體沒有關(guān)系時，就會發(fā)生這種情況磅网。如果零假設被拒絕谈截，則R≠0，這表明組內(nèi)的所有樣本對比來自不同組的任何一對樣本更相似涧偷。例如簸喂，在這種情況下，所有最相似的樣本都在樣本組內(nèi)燎潮，則R=1喻鳄。從理論上講，R<0也是可能的确封，但實際上這種情況是在生態(tài)學和微生物群研究中不太可能除呵。極端情況唉锌，R=?1，表示最相似的樣本都不在組中竿奏。第五步：檢測顯著性。與PERMANOVA檢驗一樣腥放，檢驗統(tǒng)計量R的p值是通過排列獲得的：將樣本觀察值隨機分配到組中泛啸。然后，將組內(nèi)和組間的排序相似度與隨機機會產(chǎn)生的相似度進行比較秃症。顯著性檢驗僅僅是置換后的R大于R的觀測值的分數(shù)候址，它是給定這個大小或更大的R的概率。假設檢驗背后的算法與PERMANOVA相同：如果兩組樣本單位在其物種(或其他分類群)組成上確實不同种柑，則組間的組成差異應該大于組內(nèi)的組成差異岗仑。

②基于Vegan軟件包的圖解相似度分析(ANOSIM)

相似性分析(ANOSIM)提供了一種統(tǒng)計檢驗兩組或多組抽樣單元之間是否存在顯著差異的方法。相似性分析是通過vegan包中的函數(shù)anosim()實現(xiàn)的聚请。該函數(shù)假定組內(nèi)所有排序的不同點具有大致相等的中位數(shù)和范圍荠雕。輸入數(shù)據(jù)是一個相異矩陣，可以由函數(shù)dist()或vegdist()生成驶赏。該功能還具有匯總和繪圖方法炸卑，用于執(zhí)行建模后分析。下面是一個語法示例煤傍。

anosim (data, grouping, permutations = 1000, distance = ”bray”)

其中盖文，data=數(shù)據(jù)矩陣或數(shù)據(jù)框，其中行是樣本蚯姆，列是響應變量五续，或者相異對象或相異度的對稱方陣；grouping=分組變量(因子)龄恋；permutation=用于評估ANOSIM統(tǒng)計量的重要性的排列的數(shù)目疙驾；distance=距離或相異度度量。如果輸入數(shù)據(jù)沒有不同的結(jié)構(gòu)或為對稱方陣篙挽，則需要指定距離荆萤。這里使用VDR小鼠的糞便數(shù)據(jù)來說明ANOSIM試驗。首先铣卡，我們需要加載數(shù)據(jù)和vegan包(如果它們還沒有加載)链韭。如前所述獲得分組信息。

> abund_table=read.csv("VdrFecalGenusCounts.csv",row.names=1,check.names=FALSE)

> abund_table<-t(abund_table)

> grouping<-data.frame(row.names=rownames(abund_table),t(as.data.frame (strsplit(rownames(abund_table),"_"))))

> grouping$Group <- with(grouping,ifelse(as.factor(X2)%in% c(11,12,13,14,15),c("Vdr-/-"), c("WT")))

> grouping<- grouping[,c(4)]

用Bray-Curtis相異度擬合ANOSIM：以下R碼使用Bray-Curtis相異度矩陣作為輸入數(shù)據(jù)運行ANOSIM煮落。

> library(vegan)

> bray<-vegdist(abund_table, "bray")

> anosim(bray, grouping,permutations = 1000)

Call:

anosim(dat = bray, grouping = grouping, permutations = 1000)

Dissimilarity: bray

ANOSIM statistic R: 0.19

Significance: 0.2

Permutation: free

Number of permutations: 1000

從函數(shù)ANOSIM()的輸出中可以看出敞峭，它包括建模公式以及函數(shù)調(diào)用使用的相異方法，ANOSIM統(tǒng)計量R的值蝉仇，排列的顯著性和R的排列值的個數(shù)旋讹，0.2的p值大于0.05殖蚕，表明組內(nèi)相似度不大于組間相似度，在0.05顯著水平上沉迹。我們可以得出結(jié)論睦疫，沒有證據(jù)表明組內(nèi)樣本比隨機預期的更相似。以下R碼使用大量數(shù)據(jù)框作為輸入數(shù)據(jù)運行ANOSIM鞭呕。

> anosim(abund_table, grouping, permutations = 1000, distance = "bray")

Call:

anosim(dat = abund_table, grouping = grouping, permutations = 1000,

distance = "bray")

Dissimilarity: bray

ANOSIM statistic R: 0.19

Significance: 0.19

Permutation: free

Number of permutations: 1000

ANOSIM擬合結(jié)果可以通過函數(shù)summary()進行匯總蛤育。

> fit <- anosim(bray, grouping,permutations = 1000)

> summary(fit)

Call:

anosim(dat = bray, grouping = grouping, permutations = 1000)

Dissimilarity: bray

ANOSIM statistic R: 0.19

Significance: 0.18

Permutation: free

Number of permutations: 1000

Upper quantiles of permutations (null model):

90%? 95% 97.5%? 99%

0.282 0.323 0.538 0.846

Dissimilarity ranks between and within classes:

? ? ? ? ? ? ? ? ?0% 25% 50%? 75% 100%? N

Between? ? 2? ? 9.5? ?17? ? 21.50? ?28? 15

Vdr-/-? ? ? ? 1? ? ?6.5? ?14? ? 19.75? ? 25? 10

WT? ? ? ? ? ?5? ? ? 8.0? ?11? ? ?16.50? ?22? ?3

最后，我們可以繪制結(jié)果圖(圖9.4)：

> plot(fit)

利用Jaccard相異度擬合ANOSIM：以下R碼使用Jaccard方法擬合ANOSIM葫松。

> anosim(abund_table, grouping, permutations = 1000, distance = "jaccard")

Call:

anosim(dat = abund_table, grouping = grouping, permutations = 1000,

distance= "jaccard")

Dissimilarity: jaccard

ANOSIM statistic R: 0.19

Significance: 0.2

Permutation: free

Number of permutations: 1000

用S?rensen相異度擬合ANOSIM：下列R碼適用于使用S?rensen方法的ANOSIM瓦糕。

> fit_S <- anosim(abund_table, grouping, permutations = 1000,distance ="S?rensen")

> summary(fit_S)

Call:

anosim(dat = abund_table, grouping = grouping, permutations = 1000)

Dissimilarity: bray

ANOSIM statistic R: 0.19

Significance: 0.19

Permutation: free

Number of permutations: 1000

Upper quantiles of permutations (null model):

90%? 95% 97.5%? 99%

0.282 0.344 0.546 0.846

Dissimilarity ranks between and within classes:

? ? ? ? ? ? ? ? 0%? ?25% 50%? 75% 100%? N

Between? ? 2? ? ?9.5? 17? ? ? 21.50? 28? 15

Vdr-/-? ? ? ? ?1? ? ? 6.5? 14? ? ? 19.75? 25? 10

WT? ? ? ? ? ? 5? ? ? 8.0? 11? ? ? ?16.50? 22? 3

9.4 Hypothesis Tests of Multi-response Permutation Procedures (MRPP)（多響應置換過程(MRPP)的假設檢驗）

① MRPP簡介

多響應置換過程(MRPP)是一種基于組內(nèi)和組內(nèi)差異的排列檢驗來檢驗兩組或多組樣本之間沒有差異的假設的非參數(shù)過程。測試差異可以是平均值(位置)的差異或組內(nèi)距離(擴散)的差異腋么。與PERMANOVA類似咕娄，MRPP在概念和方法上都與方差分析相結(jié)合：它比較組內(nèi)和組間的差異。潛在的想法也是一樣的珊擂。MRPP檢驗統(tǒng)計量基于組內(nèi)差異和組內(nèi)差異之間的加權(quán)平均差異圣勒。具體如下：

實施MRPP主要有五個步驟：第一步：一般情況下使用Euclidean距離計算距離矩陣，社區(qū)數(shù)據(jù)使用比例城市街區(qū)度量未玻。步驟2：計算每組Di中的平均距離灾而。步驟3：計算g組的增量。第四步：確定效果大小扳剿。統(tǒng)計量A是效應大小的度量旁趟，在ANOISIM中的解釋類似于R；它是組內(nèi)同質(zhì)性與隨機期望值的比較庇绽。

檢驗統(tǒng)計量是檢驗在零假設下不同組間是否存在差異锡搜。如果組內(nèi)異質(zhì)性等于偶然期望，則A=0瞧掺；當組內(nèi)所有項目相同時耕餐，A=1。A<0.1在生態(tài)學中普遍存在辟狈，A>0.3在生態(tài)學中較高肠缔。然而，根據(jù)我們的知識哼转，在微生物組文獻中沒有關(guān)于效應大小的報告標準明未。第五步：測試顯著性。與PERMANOVA和ANOSIM檢驗一樣壹蔓，檢驗統(tǒng)計量的p值是通過Monte Carlo排列得到的趟妥。顯著性檢驗只是簡單地計算小于觀測三角洲的排列三角洲的比例，并進行小樣本校正佣蓉。它是給定這個小的或更小的增量的概率披摄。

②以Vegan軟件包為例說明MRPP

MRPP提供了一種統(tǒng)計測試兩組或更多組采樣單元之間是否存在顯著差異的方法亲雪。MRPP是通過vegan包中的函數(shù)mrpp()實現(xiàn)的。如果輸入數(shù)據(jù)不同疚膊，則可以直接使用的义辕。如果它是具有按響應進行觀測的數(shù)據(jù)結(jié)構(gòu)的矩陣，則在使用之前需要通過函數(shù)vegdist()計算相異度寓盗。默認距離是歐幾里得方法终息，但是還應用了vegdist()中的其他不同之處。該函數(shù)具有匯總和繪圖兩種方法贞让，用于執(zhí)行建模后分析。以下是一個示例語法：mrpp(data, grouping, permutations = 1000, distance = “bray”)

其中柳譬，data=數(shù)據(jù)矩陣或數(shù)據(jù)框喳张，其中行是觀測值，列是響應美澳，或者不同的對象或不同的對稱方陣销部。響應可以是單變量或多變量舅桩；grouping=分組變量(因子)雨膨；排列=用于評估MRPP統(tǒng)計量重要性的排列數(shù)目擂涛；距離=距離或相異度量撒妈。如果輸入數(shù)據(jù)沒有不同的結(jié)構(gòu)或為對稱方陣，則需要指定距離排监。與ANOSIM測試一樣，我們需要訪問數(shù)據(jù)并加載vegan包棋蚌。

用Bray-Curtis相異度擬合MRPP：下面的R碼使用Bray-Curtis相異度矩陣作為輸入數(shù)據(jù)來運行MRPP谷暮。

> mrpp(bray, grouping,permutations = 1000)

Call:

mrpp(dat = bray, grouping = grouping, permutations = 1000)

Dissimilarity index: bray

Weights for groups:? n

Class means and counts:

? ? ? ? ?Vdr-/- WT?

delta 0.44? 0.427

n? ? ? ? ?5? ? ? 3? ?

Chance corrected within-group agreement A: 0.0502

Based on observed delta 0.4353 and expected delta 0.4583

Significance of delta: 0.17

Permutation: free

Number of permutations: 1000

觀測delta和預期delta分別為0.4353和0.4583瞒瘸。差值的顯著性為0.0502.13情臭，校正后的組內(nèi)協(xié)議A的機率為0.。在樣本量較小的情況下娃惯，我們得出的結(jié)論是趾浅，在0.05的顯著性水平上皿哨，兩個屬群之間沒有統(tǒng)計上的顯著差異纽谒。以下代碼使用豐度數(shù)據(jù)框作為輸入數(shù)據(jù)運行MRPP：

> mrpp(abund_table, grouping, permutations = 1000, distance = "bray")

Call:

mrpp(dat = abund_table, grouping = grouping, permutations = 1000,distance= "bray")

Dissimilarity index: bray

Weights for groups:? n

Class means and counts:

? ? ? ? Vdr-/- WT?

delta 0.44? 0.427

n? ? ? ? 5? ? ? 3? ?

Chance corrected within-group agreement A: 0.0502

Based on observed delta 0.4353 and expected delta 0.4583

Significance of delta: 0.15

Permutation: free

Number of permutations: 1000

使用Meandist()獲取平均距離矩陣：函數(shù)meandist()計算平均簇內(nèi)(塊)相異度(對角線)和簇間(塊)相異度(非對角線元素)的矩陣证膨，以及分組計數(shù)的屬性n央勒。

> meandist(bray, grouping,permutations = 1000)

? ? ? ? Vdr-/-? WT

Vdr-/- 0.4401 0.4766

WT? ? 0.4766 0.4273

attr(,"class")

[1] "meandist" "matrix"?

attr(,"n")

grouping

Vdr-/-? WT

5? ? ? 3

使用meandist()和Summary()獲取平均距離和匯總統(tǒng)計信息：我們可以對函數(shù)meandist()對象使用函數(shù)Summary()來查找這些不同之處的類內(nèi)崔步、類間和總體平均值井濒，以及具有所有權(quán)重類型選項和分類強度的MRPP統(tǒng)計數(shù)據(jù)眼虱。

> bray_mrpp <- meandist(bray, grouping,permutations = 1000,distance ="bray",weight.type = 1)

> summary(bray_mrpp)

Mean distances:

? ? ? ? ? ? ? ? ? ? ? ? ? Average

within groups? ? ? 0.4372

between groups? 0.4766

overall? ? ? ? ? ? ? ? 0.4583

Summary statistics:

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Statistic

MRPP A weights n? ? ? ? ? 0.05016

MRPP A weights n-1? ? ? ? 0.04900

MRPP A weights n(n-1)? ? 0.04612

Classification strength? 0.04131

9.5 Compare Microbiome Communities Using the GUniFrac Package（使用GUniFrac軟件包比較微生物群落）

① UniFrac、加權(quán)UniFrac和廣義UniFrac距離度量簡介

在多變量微生物群落組成分析中过牙，未加權(quán)和加權(quán)UniFrac距離的度量方法被最廣泛地應用于兩個微生物群落的度量寇钉。UniFrac距離扫倡，也稱為未加權(quán)的UniFrac距離撵溃，是由Lozupone等人提出的缘挑。這是一種基于系統(tǒng)發(fā)育距離來估計微生物組樣本之間差異的方法语淘。UniFrac距離度量的目標是實現(xiàn)來自不同條件的微生物樣本之間的客觀比較惶翻。未加權(quán)的UniFrac定義如下：

按照這種定義畔咧，UniFrac度量的計算方法是將兩個樣本之間沒有共享的分支長度除以任一樣本覆蓋的分支長度，而不是兩個樣本都覆蓋的分支長度。距離0表示兩個樣本相同垃沦，距離1表示兩個樣本沒有共同的分類群肢簿。作為缺失的二元檢驗，未加權(quán)的UniFrac距離僅考慮分類單元的存在和缺失信息收夸；它在檢測稀有譜系的豐度變化方面是最有效的卧惜。然而咽瓷，這種距離度量的定義完全忽略了類群豐度信息茅姜。在原始的未加權(quán)方法中增加了比例加權(quán)钻洒，因此將這種新的UniFrac度量稱為加權(quán)UniFrac航唆。加權(quán)UniFrac距離利用類群豐度信息糯钙，并根據(jù)豐度差對分支長度進行加權(quán)任岸；在加權(quán)UniFrac距離度量中享潜，系統(tǒng)發(fā)育樹的每個分支長度是根據(jù)兩個樣本之間分類單元豐度比例的差異來加權(quán)的剑按，而不是只看是否有分類單元艺蝴。加權(quán)UniFrac距離被定義為：

該定義的表述有幾個特征：(1)它用豐度差來加權(quán)分支長度姑荷；(2)即使我們將豐度數(shù)據(jù)轉(zhuǎn)換成有無數(shù)據(jù)鼠冕，dW也不能被簡化為dU懈费；以及(3)dW使用絕對比例差楞捂，這導致DWI的值主要由比例大的分支決定寨闹，并且對小比例的分支上的豐度變化不那么敏感繁堡。通過將比例加權(quán)添加到UniFrac距離椭蹄，加權(quán)的UniFrac距離減少了低豐度分類群被表示為0或由取決于采樣深度的低計數(shù)表示的問題罩润。在加權(quán)的UniFrac中割以，低豐度分類群的重量要低得多严沥，因此對度量報告的總距離的影響較小消玄。因此翩瓜，加權(quán)的UniFrac既可以檢測每個譜系中存在多少序列的變化兔跌，也可以檢測存在哪些分類群的變化；檢測豐富血統(tǒng)的變化是最敏感的层亿。然而，無論是未加權(quán)的還是加權(quán)的UniFrac距離在檢測適度豐富的譜系中的變化方面都可能不是非常強大碌更，因為他們太看重稀有血統(tǒng)或最豐富的血統(tǒng)痛单。因此提出了以下廣義UniFrac距離來統(tǒng)一加權(quán)UniFrac距離和未加權(quán)UniFrac距離旭绒。

② Breast Milk Data Set（母乳數(shù)據(jù)集）

母乳數(shù)據(jù)集來自最近發(fā)表的兩項研究重父。它收集了58個微生物組樣本房午，這些樣本取自正在哺乳的加拿大高加索女性郭厌。母乳是發(fā)育中的嬰兒的重要細菌來源，會影響新生兒的細菌組成液走，進而影響日后的疾病風險缘眶。在第一篇論文中巷懈，作者使用這個數(shù)據(jù)集來比較早產(chǎn)和足月產(chǎn)顶燕、剖腹產(chǎn)(選擇性和非選擇性)和陰道分娩以及男嬰和女嬰之間的細菌圖譜。在第二期出版物中恳谎，使用相同的數(shù)據(jù)集來說明作者自己開發(fā)的工具因痛，包括未加權(quán)的UniFrac鸵膏、加權(quán)的UniFrac较性、informationUniFrac和Ratio? UniFrac责循。在這里院仿，我們使用這個數(shù)據(jù)集來比較使用GUniFrac軟件包的微生物群落歹垫。

③?使用GUniFrac軟件包比較微生物群落

軟件包GUniFrac是用來實現(xiàn)廣義的uniFrac距離來比較微生物群落的排惨。作為PERMANOVA過程，在GUniFrac軟件包中辕宏，基于排列來評估偽F統(tǒng)計量的重要性瑞筐。如果一種環(huán)境所特有的樹木的比例比偶然預期的要大聚假，那么這兩個群落就被認為是不同的。在此軟件包中闯袒，還可以實現(xiàn)未加權(quán)和加權(quán)UniFrac以及方差調(diào)整后的加權(quán)UniFrac距離。

其中mi是第i個分支上來自兩個社區(qū)的個人/讀取總數(shù)，m是個人/讀取總數(shù)访锻。它的開發(fā)是為了說明這樣一個事實期犬，即加權(quán)UniFrac距離沒有考慮隨機抽樣下權(quán)重的變化璃谨，從而導致檢測群落之間差異的能力較小(Chang等人佳吞。GUniFrac包依賴于vegan包和ape包底扳，作者還建議使用包ade4。下面給出一種用法算芯。 GUniFrac(otu.tab, tree, alpha = c(0, 0.5, 1))

其中熙揍，out.tab=具有行(=樣本)和列(=OTU)的OTU計數(shù)表届囚；tree= rooted phylogenetic tree of R class “phylo”；alpha= parameter controlling weight on abundant lineages蛔添。UniFrac測量計算需要兩條信息：計數(shù)表和系統(tǒng)發(fā)育樹迎瞧。因此，運行GUniFrac軟件包需要兩個數(shù)據(jù)集：計數(shù)表(如OTU表)和系統(tǒng)發(fā)生樹足绅。在下面氢妈，我們將一步一步地說明如何使用GUniFrac包來比較微生物群落厕怜。