久久国产成人av_抖音国产毛片_a片网站免费观看_A片无码播放手机在线观看,色五月在线观看,亚洲精品m在线观看,女人自慰的免费网址,悠悠在线观看精品视频,一级日本片免费的,亚洲精品久,国产精品成人久久久久久久

分享

三萬字懂得微生物群落差異分析(下)

 微生信生物 2021-03-12

9.3  用ANOSIM假設(shè)檢驗(yàn)組間的不同

9.3.1相似性分析介紹(ANOSIM)

ANOSIM僅僅是基于兩個(gè)距離矩陣之間標(biāo)準(zhǔn)化秩相關(guān)的Mantel測試的修正版本,。

ANOSIM檢驗(yàn)是由Clarke(1993)開發(fā)的一種無分布的多元數(shù)據(jù)分析方法,經(jīng)常被群落生態(tài)學(xué)家和微生物組研究人員使用,。它是一種非參數(shù)程序,用于檢驗(yàn)兩組或多組樣本之間無差異的假設(shè),,基于組間和組內(nèi)相似性的排列檢驗(yàn)(Clarke 1993),。它比較了物種(或任何其他分類群)豐度和組成在采樣單元(如:就分組因素或?qū)嶒?yàn)處理水平而言,與ANOVA分析類似,,ANOSIM將組成員資格或處理水平作為因素,,并將其建模作為解釋變量。相似度的分析基于一個(gè)簡單的想法:如果測試的組是有意義的,,那么組內(nèi)的樣本應(yīng)該比不同組的樣本在組成上更相似,。因此,零假設(shè)是:在測試組或治療條件的成員之間沒有差異,。本實(shí)驗(yàn)采用Bray-Curtis相似性測度法,。

anosim檢驗(yàn)統(tǒng)計(jì)量是基于組間和組內(nèi)的平均秩差。如下:

里面R是檢驗(yàn)統(tǒng)計(jì)量,,是相對群內(nèi)差異的一個(gè)指標(biāo),。

M=N(N-1)/2  是樣品對數(shù)。

N是為樣本總數(shù)(受試者),。

rB是組間相似性排序的平均值

rW是組內(nèi)相似性排序的平均值,。

執(zhí)行ANOSIM有五個(gè)主要步驟:

步驟1:計(jì)算距離矩陣。

步驟2:計(jì)算秩差值,,將最小的不相似值賦1,。

步驟三:計(jì)算組間和組內(nèi)秩差的平均值。

步驟4:用上式計(jì)算檢驗(yàn)統(tǒng)計(jì)量R

R被解釋為一個(gè)相關(guān)系數(shù),,作為其他類型的相關(guān)系數(shù),,如Pearson’s相關(guān)系數(shù),是“效應(yīng)大小”的衡量標(biāo)準(zhǔn),。檢驗(yàn)統(tǒng)計(jì)量是為了檢驗(yàn)在原假設(shè)下組間是否存在差異,。如果零假設(shè)是正確的,那么R=0,,這意味著平均而言,,組間和組內(nèi)的差異是相同的。當(dāng)高相似性和低相似性完全混合且與群體沒有關(guān)系時(shí),就會發(fā)生這種情況,。如果零假設(shè)被拒絕,,那么R不等于0,,這表明組內(nèi)所有的樣本對比任何來自不同組的樣本對更相似,。例如,在這種情況下,,所有最相似的樣本都在樣本組內(nèi),,則R = 1。理論上,,R < 0也是有可能的,。但實(shí)際上這種情況是在生態(tài)和微生物組研究中不太可能。極端情況R =-1,,表示最相似的樣本都不在組內(nèi),。

步驟五4:檢驗(yàn)顯著性。

正如在PERMANOVA檢驗(yàn)中,,檢驗(yàn)統(tǒng)計(jì)量R的p值是通過排列獲得:將樣本觀察值隨機(jī)分配給組,。然后,將組內(nèi)和組間的排序相似度與隨機(jī)生成的相似度進(jìn)行比較,。顯著性檢驗(yàn)簡單地說就是R大于R的觀測值的排列R的分?jǐn)?shù),。它是給定R大于或等于這個(gè)值的概率。假設(shè)檢驗(yàn)背后的算法與PERMANOVA相同:如果兩組抽樣單位在其物種(或其他分類群)組成上確實(shí)不同,,那么組間的組成差異應(yīng)該大于組內(nèi)的組成差異,。

9.3.2  用vegan包說明ANOSIM的相似性分析

相似性分析(ANOSIM)提供了一種統(tǒng)計(jì)檢驗(yàn)兩組或多組抽樣單位之間是否存在顯著差異的方法。相似性分析是通過函數(shù)anosim()在vegan包中實(shí)現(xiàn)的,。該函數(shù)假設(shè)所有組內(nèi)的等級差異的中位數(shù)和范圍都是相等的,。輸入數(shù)據(jù)是一個(gè)不同矩陣。它可以由函數(shù)dist()或vegdist()生成,。該功能還具有總結(jié)和繪圖方法來進(jìn)行后期建模分析,。下面是一個(gè)語法示例。

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

其中,,data為數(shù)據(jù)矩陣,,其中行為樣本,列為響應(yīng)變量,?;虿煌奈矬w或?qū)ΨQ的不同方陣;; grouping = grouping variable (一個(gè)因素); permutations = 評估ANOSIM統(tǒng)計(jì)量的顯著性的排列數(shù); distance = 距離或不相似度量。如果輸入數(shù)據(jù)不具有不相似結(jié)構(gòu)或?yàn)閷ΨQ方陣,,則需要指定距離,。這里用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的不同點(diǎn)來匹配ANOSIM

以下R代碼使用Bray-Curtis不同矩陣作為輸入數(shù)據(jù)運(yùn)行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)用中使用的不同方法。由排列值和排列值個(gè)數(shù)的顯著性可知,,0.2的p值大于0.05,,表明在0.05顯著水平下,組內(nèi)相似度不大于組間相似度,。我們可以得出結(jié)論,,沒有證據(jù)表明組內(nèi)樣本比隨機(jī)概率預(yù)期的更相似

下面的R代碼使用豐度數(shù)據(jù)幀作為輸入數(shù)據(jù)運(yùn)行ANOSIM。

> anosim(abund_table, grouping, permutations = 1000, distance = "bray")
distance = "bray")
Dissimilarity: bray
ANOSIM statistic R: 0.19
Significance: 0.19
Permutation: free
Number of permutations: 1000

ANOSIM的匹配結(jié)果可以由summary().函數(shù)概括,。

> 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. 以下使用S?rensen方法適用于ANOSIM rensen

> 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  Multi-response Permutation Procedures的假設(shè)檢驗(yàn)

9.4.1 介紹MRPP

多響應(yīng)置換程序(MRPP)是一種非參數(shù)程序,,用于檢驗(yàn)兩組或多組樣本之間無差異的假設(shè),基于組間和組內(nèi)差異的置換檢驗(yàn)(Mielke 1984,,1991),。測試差異可能會使均值(位置)或組內(nèi)距離差異(spread) (Warton et al. 2012)。

與PERMANOVA相似,,在概念和方法上,,MRPP與方差分析(ANOVA)相結(jié)合:它比較組內(nèi)和組間的差異。潛在的原理也是一樣的,。MRPP檢驗(yàn)統(tǒng)計(jì)量是基于組間和組內(nèi)差異的加權(quán)均值之差,。它如下所示:

式中,delta= MRPP檢驗(yàn)統(tǒng)計(jì)量,,為抽樣單位間兩兩不相似的組內(nèi)均數(shù)的總體加權(quán)均值:Ci=ni/N,;N=項(xiàng)目總數(shù): n;=組i中項(xiàng)目的數(shù)量。

實(shí)施MRPP主要有五個(gè)步驟: 第一步:一般情況下使用歐幾里得距離計(jì)算距離矩陣 第二步:計(jì)算各組di的平均距離,。第三步:計(jì)算g多組的增量 第四步:確定效果大小,。統(tǒng)計(jì)量A是對效應(yīng)大小的衡量,它被解釋為ANOISIM中的R;它是作為組內(nèi)同質(zhì)性與隨機(jī)期望的比較而得到的,。

檢驗(yàn)統(tǒng)計(jì)量是檢驗(yàn)在零假設(shè)下組間無差異,。如果組內(nèi)異質(zhì)性隨機(jī)等于期望,則A=0:當(dāng)組內(nèi)所有項(xiàng)各都相同時(shí),,則A= 1,。A < 0.1在生態(tài)學(xué)中很常見;> 0.3在生態(tài)學(xué)中是相當(dāng)高的。然而,,據(jù)我們所知,,目前尚無微生物組效應(yīng)大小的報(bào)道標(biāo)準(zhǔn),。

第五步:檢驗(yàn)顯著性。正如在PERMANOVA和ANOSIM檢驗(yàn)中,,檢驗(yàn)統(tǒng)計(jì)量delta的p值是通過Monte Carlo置換得到的,。顯著性檢驗(yàn)是簡單地用小于觀測到的增量的置換增量的比例,并進(jìn)行小樣本校正,。它是在這個(gè)或更小的條件下的概率,。

9.4.2 用vegan包做MRPP

MRPP提供了一種統(tǒng)計(jì)檢驗(yàn)兩組或兩組以上采樣單位之間是否存在顯著差異的方法。MRPP是通過vegan包中的mrpp()函數(shù)實(shí)現(xiàn)的,。如果輸入數(shù)據(jù)不一致,,則可以直接使用,。如果它是一個(gè)由響應(yīng)數(shù)據(jù)結(jié)構(gòu)觀察到的矩陣,,那么在使用之前需要由函數(shù)vegdist()計(jì)算不相似性。默認(rèn)的距離是歐幾里得方法,,但其他不同的vegdist()也被應(yīng)用,。該函數(shù)有總結(jié)和繪圖方法來進(jìn)行后期建模分析。下面是一個(gè)語法示例:

mrpp(data, grouping, permutations = 1000, distance = “bray”)

其中,,data =數(shù)據(jù)矩陣或數(shù)據(jù)框架,,其中行是觀測值,列是結(jié)果,,或一個(gè)不相似對象或?qū)ΨQ的不相似方陣,。結(jié)果可能是單變量的或多變量的;grouping =分組變量(一個(gè)因素);permutations =評估MRPP統(tǒng)計(jì)量的顯著性的置換數(shù):distance = distance or dissimilarity measure如果輸入數(shù)據(jù)沒有dissimilarity structure或者是一個(gè)對稱的方陣,那么需要指定距離,。在ANOSIM測試中,,我們需要訪問數(shù)據(jù)并加載vegan包。

用Bray-Curtis的不同點(diǎn)來匹配MRPP以下R代碼使用Bray-Curtis不同矩陣作為輸入數(shù)據(jù)運(yùn)行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

觀測到的和預(yù)期的deltas分別為0.4353和0.4583,。delta的顯著性為0.13,,組內(nèi)一致性A校正的概率為0.0502。由于樣本量較小,,我們認(rèn)為在0.05的顯著水平下,,兩個(gè)屬之間的差異不具有統(tǒng)計(jì)學(xué)意義。

以下代碼使用豐度數(shù)據(jù)幀作為輸入數(shù)據(jù)運(yùn)行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ù)獲得平均距離矩陣函數(shù)meandist()計(jì)算群內(nèi)(塊)差異(對角線)和簇間(塊)差異(非對角線元素)的均值矩陣,,以及分組計(jì)數(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()函數(shù)獲得平均距離和匯總統(tǒng)計(jì)我們可以使用meandist()summary()函數(shù)來查找組內(nèi)、組間和總體方法以及他們的差異,,以及所有重要的MRPP統(tǒng)計(jì)數(shù),。選項(xiàng)和分類強(qiáng)度。

> 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使用GUniFrac包比較微生物群落

9.5.1  介紹UniFrac,,加權(quán)UniFrac和廣義UniFrac距離矩陣

在多變量微生物組組成分析中,,加權(quán)和未加權(quán)的 UniFrac distance是測量兩種微生物組群體的的最常用方法,。UniFrac距離,也被稱為未加權(quán)UniFrac距離,,由Lozupone等人在2005年提出(Lozupone and Knight 2005),。它是根據(jù)系統(tǒng)發(fā)育距離估計(jì)微生物組樣品之間差異的一種方法。UniFrac距離度量的目標(biāo)是使來自不同條件的微生物組樣本之間的客觀比較成為可能,。unweighted UniFrac定義如下:

式中,,d = unweighted UniFrac distance;A,B=微生物組群落A與B;n=系統(tǒng)發(fā)育樹的根狀分支;bi=i枝長度,piA和piB=群落A和B從分支i下降的類群比例。I()是指示物種是否存在于i支的指示函數(shù),。在上式中,,群落A和群落B從i支下降的類群比例分別定義為>0和I(piA>0)和I(piB>0)。

按照這種定義,,UniFrac度量是通過將兩個(gè)樣本之間不共享的分支長度除以其中一個(gè)樣本(但不是兩個(gè)樣本都包含的分支長度)來計(jì)算的(Lozupone和Knight 2005),。距離為0表示兩個(gè)樣本完全相同,距離為1表示兩個(gè)樣本沒有共同類群,。unweighted UniFrac distance作為缺席二值檢驗(yàn),,只考慮類群的存在和缺席信息;它是檢測稀有血統(tǒng)豐度變化的最有效方法。然而,,這種距離度量的定義完全忽略了類群豐度信息(Chen et al .2012).2007年,,Lozupone等人在原始的未加權(quán)方法(Lozupone et al. 2007)的基礎(chǔ)上增加了比例加權(quán),因此將這種新的UniFrac方法稱為加權(quán)UniFrac,。加權(quán)UniFrac距離利用類群豐度信息,,并對豐度差異的分支長度加權(quán);在加權(quán)的UniFrad距離度量中,系統(tǒng)發(fā)育樹的每個(gè)分支長度的加權(quán)是由兩個(gè)樣本之間的類群豐度比例的差異,,而不是只看類群的存在或不存在,。加權(quán)UniFrac距離(Lozupone et al. 2007)定義為

式中,dW=(標(biāo)準(zhǔn)化)加權(quán)UniFrac距離;A,、B=微生物組群落A,、B;n =系統(tǒng)發(fā)育樹的根狀分支;bi =分支i的長度。

該定義的公式有以下幾個(gè)特點(diǎn):(1)以豐度差異為分支長度加權(quán);==(2)即使將豐度數(shù)據(jù)轉(zhuǎn)化為存在缺失數(shù)據(jù),,也不能將dW簡化為du;(3) dw采用絕對比例差,,導(dǎo)致dw值主要由大比例分支決定,對小比例分支豐度變化不太敏感(Chen et al. 2012),。通過向UniFrac距離添加一個(gè)比例加權(quán),,加權(quán)的UniFrac距離減少了低豐度分類群被表示為0或根據(jù)采樣深度的低計(jì)數(shù)的問題。在加權(quán)的UniFrac中,,低豐度分類群的權(quán)重更低,,因此對度量報(bào)告的總距離的影響更小(Wong et al. 2016)。因此,,加權(quán)的UniFrac可以檢測每個(gè)譜系中存在的序列數(shù)量和分類群的變化(Lozupone et al. 2007);它對檢測豐富基因的變化最為敏感(Chen et al. 2012),。然而,,未加權(quán)或加權(quán)的UniFrac距離在檢測中等豐富譜系中的變化時(shí)可能都不是很有效(Chen et al. 2012),因?yàn)樗鼈儗币娮V系或最豐富的譜系都賦予了過多的權(quán)重,。因此,,Chen等人提出了以下廣義UniFrac距離,以統(tǒng)一加權(quán)的UniFrac距離和非加權(quán)的UniFrac距離,。

式中,d(α)=廣義UniFrac距離,;α∈[0,1] 是用來控制來自高豐度分支的貢獻(xiàn);A,B=分別是微生物群落A,B;n=生根的系統(tǒng)發(fā)育樹的分支;bi=分支i的長度,;∑nbi(piA+piB)(α)=標(biāo)準(zhǔn)因子以確保d(α)值在[0,1].廣義UniFrac 距離包含一個(gè)額外的參數(shù)α控制家系豐度的權(quán)重,,所以距離不由高豐度家系決定。α=0.5已經(jīng)遠(yuǎn)超于最大值,。UniFrac距離的標(biāo)準(zhǔn)化公式有幾個(gè)特征:(1)利用相對差值|piA-piB|/(piA+piB)對大比例分支的權(quán)重進(jìn)行了衰減,,權(quán)重值在[0,1]區(qū)間(2)根據(jù)這個(gè)公式,如果我們將豐度數(shù)據(jù)轉(zhuǎn)換為存在/缺失數(shù)據(jù),,現(xiàn)在dW可以簡化為dU,,當(dāng)我們令式(9.11)中α=1,,d(α)減少到dW,。

當(dāng)α=0,我們得到

基于UniFrac距離的偽f統(tǒng)計(jì)量定義為

其中tr(.)是矩陣的跡函數(shù),,H=X(XTX)-1XT  為設(shè)計(jì)矩陣X的hat(投影)矩陣,,G為高爾居中矩陣,n為樣本個(gè)數(shù),,m為預(yù)測因子個(gè)數(shù),。設(shè)dij為群落i與j之間的廣義UniFrac距離,且表示為A=(aij)=(-1/2dij2).高爾矩陣定義為G=(I-11/n)A(I-11/n),。未加權(quán),、加權(quán) UniFrac距離和廣義UniFrac距離的詳細(xì)的檢驗(yàn)統(tǒng)計(jì)量及統(tǒng)一公式,讀者可以參考作者的原始論文,。

9.5.2  母乳數(shù)據(jù)集

母乳數(shù)據(jù)集來自兩項(xiàng)最近發(fā)表的研究(Urbaniak et al.),。2016;Wong等人2016年)。它從哺乳期收集了58個(gè)微生物樣本加拿大白人女性,。母乳是發(fā)育中的嬰兒的重要細(xì)菌來源,,并影響新生兒的細(xì)菌組成,進(jìn)而影響到以后生活中的疾病風(fēng)險(xiǎn),。在第一篇論文中,,作者使用這組數(shù)據(jù)比較了早產(chǎn)和出生的、剖腹產(chǎn)(選擇性和非選擇性),、陰道分娩以及男嬰和女嬰之間的細(xì)菌分布,。在第二篇文章中,,使用相同的數(shù)據(jù)集來說明作者自己開發(fā)的工具,包括unweighted UniFrac,、weighted UniFrac,、information UniFrac和ratio UniFrac。在這里,,我們使用這個(gè)數(shù)據(jù)集來比較GUniFrac軟件包中的微生物群落,。

9.5.3

使用GUniFrac包比較微生物群落

GUniFrac包被發(fā)展以完成廣義的uniFrac距離來比較生物群落,正如PERMANOVA程序,, GUniFrac包基于置換估計(jì)偽F的統(tǒng)計(jì)值的顯著性,,如果樹在一個(gè)環(huán)境中所占的比例大于偶然的預(yù)期值,那么這兩個(gè)群落就被認(rèn)為是不同的,,在這個(gè)包中,,加權(quán)和不加權(quán)UniFrac,以及方差調(diào)整加權(quán)UniFrac距離也可以實(shí)現(xiàn),, UniFrac距離可以被計(jì)算出來,,并且在GUniFrac包中比較微生物群落,包括dW,d(0.5),d(0),d(U),以及方差以調(diào)整加權(quán)UniFrac距離dVAW,,以上部分,我們介紹了dW,d(0)和d(U)的公式,,d(0.5)代表一系列距離的中間值,通過以下公式計(jì)算出來:

方差調(diào)整加權(quán)UniFrac距離dVAW計(jì)算公式:

其中mi為第i支兩個(gè)社團(tuán)的總個(gè)體數(shù)/reads, m為總個(gè)體數(shù)/reads,。它的開發(fā)是為了考慮這樣一個(gè)事實(shí),,即加權(quán)UniFrac距離不考慮隨機(jī)抽樣下權(quán)重的變化,從而導(dǎo)致較少的功率檢測社區(qū)之間的差異(Chang et al. 2011),。GUniFrac包依賴于vegan和ape包,,作者還建議使用ade4包。一個(gè)用法如下所示:

GUniFrac(otu.tab, tree, alpha = c(0, 0.5, 1))

UniFrac測量計(jì)算需要兩條信息:嘔吐表和系統(tǒng)發(fā)育樹,。因此,,運(yùn)行g(shù)uunifrac包需要兩個(gè)數(shù)據(jù)集:一個(gè)計(jì)數(shù)表,如OTU表和系統(tǒng)發(fā)育樹,。接下來,,我們將逐步演示如何使用GUniFrac軟件包來比較微生物群落。

第一步:加載并讀取合適格式的OTU表設(shè)置目錄的方法與我們在前面章節(jié)中所做的相同,?!皌d_OTU_tag_mapped_lineage.txt”數(shù)據(jù)集包括兩條信息:最后一列的OTU計(jì)數(shù)表和分類。來自原始數(shù)據(jù)集作者的部分R代碼是有用的,。我們在這里修改它們,,并解釋和評論代碼的作用。

> setwd("E:/Home/MicrobiomeStatUsingR/Analysis/")
> otu_tab <- read.table("td_OTU_tag_mapped_lineage.
txt", header=T, sep="\t", row.names=1, comment.
char="", check.names=FALSE)

在參數(shù)中,,row.names = 1 參數(shù)指定表中包含行名的列1;comment.char  = “ “用于完全關(guān)閉注釋的解釋,。它比默認(rèn)的read.table()快得多,。數(shù)據(jù)集的前6行“td_OTU_tag_mapped_lineage.txt”的內(nèi)容如下:

> head(otu_tab)
S31 S1 S42 S13_T2 S30 S50 S43 S20 S29 S47U S26 S13_T3 S33 S8L
0 38 36 30 14 13 18 27 38 49 7 5251 12 17 35
1 2866 15069 42985 3292 1223 1056 3959 3021 7023 1856 7993 2571 13306 1101
2 437 9831 4628 4231 4473 718 3843 441 4311 2496 5672 4940 4112 4830
3 3356 7407 5355 62 121 14 3616 6184 334 108 31 30 13 56
4 12 478 16 342 23 8 18 97 390 126 4 271 53 7
5 145 238 234 3 56 1 109 156 374 28 0 45 0 1
……
taxonomy
0
Bacteria;Proteobacteria;Gammaproteobacteria;Pasteurellales;Pasteurellaceae;Pasteurella;|93
1
Bacteria;Firmicutes;Bacilli;Bacillales;Staphylococcaceae;Staphylococcus;|77
2
Bacteria;Proteobacteria;Gammaproteobacteria;Pseudomonadales;Pseudomonadaceae;Pseudomonas;|92
3 Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacteriales;Enterobacteriaceae;EscherichiaShigella;|98
4
Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacteriales;Enterobacteriaceae;Klebsiella;|75
5
Bacteria;Actinobacteria;Actinobacteria;Corynebacteriales;Corynebacteriaceae;Corynebacterium;|88

然而,為了正確使用guunifrac包,,OTU表應(yīng)該是一個(gè)數(shù)字矩陣,。下面的R代碼用于從原始數(shù)據(jù)集中刪除分類列。

> taxonomy <- otu_tab$taxonomy
> otu_tab <- otu_tab[-length(colnames(otu_tab))]

otu表需要通過otu矩陣轉(zhuǎn)置為樣本,。

> otu_tab <- t(as.matrix(otu_tab))
> head(otu_tab)
......
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
S31 38 2866 437 3356 12 145 43 0 20 19 33 27 16 21 79 40 30 12
S1 36 15069 9831 7407 478 238 1368 25 306 2821 91 169 323 475 1290 3795 532 186
S42 30 42985 4628 5355 16 234 14491 2 357 153 48 109 136 520 781 617 727 867
S13_T2 14 3292 4231 62 342 3 282 1 188 51 64 168 161 574 557 33 408 7
S30 13 1223 4473 121 23 56 3 0 123 3 13 3 279 207 703 123 617 5
S50 18 1056 718 14 8 1 3 1 20 5 0 1 35 49 119 28 83 5

第二步:使用抽平將樣品OTU計(jì)數(shù)標(biāo)準(zhǔn)化到一個(gè)標(biāo)準(zhǔn)測序深度UniFrac的作者推薦使用稀疏數(shù)據(jù)(Carcer等人,,2011年),并將其用于計(jì)算未加權(quán)的UniFrac距離矩陣,,而非稀疏數(shù)據(jù)則使用自定義的UniFrac腳本用于加權(quán),、信息和比例UniFrac (Wong等人,2016年),。然而,,在GUuniFrac軟件包手冊中,所有的UniFrac距離計(jì)算都使用了稀有數(shù)據(jù),。因?yàn)镚UuniFrac對不同的測序深度很敏感,,為了在相同的基礎(chǔ)上比較微生物群,可以使用稀疏法(Chen 2012),。因此,,我們遵循包手冊;在執(zhí)行未加權(quán),加權(quán)UniFrac和方差調(diào)整加權(quán)UniFrac距離前細(xì)化樣本 ,??梢允褂胿egan package中的rrarefy()或rarefy()函數(shù)來實(shí)現(xiàn)。

> library(vegan)
> otu_tab_rarefy <- rrarefy(otu_tab, min(apply(otu_tab,1,sum)))

步驟3:讀取系統(tǒng)發(fā)育樹GUniFrac包需要一個(gè)根樹作為輸入數(shù)據(jù),。我們可以使用phangorn包中的函數(shù)midpoint()來獲取根樹。

> library(phangorn)
> otu_tree <- midpoint(otu_tree)
> otu_tree
Phylogenetic tree with 115 tips and 114 internal nodes.
Tip labels:
1, 11, 6686, 18, 230, 82, ...
Node labels:
NA, 71.9, 79.1, 86.9, 87.1, 52.4, ...
Rooted; includes branch lengths.

第四步:計(jì)算UniFracs現(xiàn)在,,可以使用GUuniFrac軟件包計(jì)算UniFracs

> library(GUniFrac)
> #Calculate the UniFracs
> unifracs <- GUniFrac(otu_tab_rarefy, otu_tree, alpha=c(0, 0.5, 1))
$unifracs
> dw <- unifracs[,, "d_1"] # Weighted UniFrac
> du <- unifracs[,, "d_UW"] # Unweighted UniFrac
> dv <- unifracs[,, "d_VAW"]# Variance adjusted weighted UniFrac
> d0 <- unifracs[,, "d_0"] # GUniFrac with alpha 0
> d5 <- unifracs[,, "d_0.5"]# GUniFrac with alpha 0.5

第五步:進(jìn)行PERMANOVA比較UniFrac的一項(xiàng)測量為了檢驗(yàn)UniFrac測量的假設(shè),,需要從元數(shù)據(jù)中提取群組信息。我們使用read.table()函數(shù)讀取元數(shù)據(jù),,如下所示:

> meta_tab<- read.table("metadata.txt", header=T, sep="\t",
row.names=1, comment.char="", check.names=FALSE)

匹配OTU表和分組數(shù)據(jù),,以保留只出現(xiàn)在兩個(gè)數(shù)據(jù)集中的樣本.

> otu_meta_matched <- match(rownames(meta_tab),rownames(otu_tab))
> otu_meta_matched <- otu_meta_matched[!is.na(otu_meta_matched)]
> otu_tab <- otu_tab[otu_meta_matched,]
> meta_tab_Ordered <- meta_tab[match(rownames(otu_tab),rownames
(meta_tab)),]

下面的R代碼使用了vegan package中的adonis()函數(shù)來執(zhí)行 PERMANOVA比較妊娠距離系列的中間距離。感興趣的讀者可以嘗試UniFrac的其他距離指標(biāo),。

> set.seed(123)
> adonis(as.dist(d5) ~ meta_tab$Gestation)
Call:
adonis(formula = as.dist(d5) ~ meta_tab$Gestation)
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
meta_tab$Gestation 4 0.4336 0.108391 1.2001 0.08305 0.22
Residuals 53 4.7869 0.090319 0.91695
Total 57 5.2205 1.00000

類似地,,性別和奶類型的影響可以被測試如下。

> adonis(as.dist(d5) ~ meta_tab$Gender)
Call:
adonis(formula = as.dist(d5) ~ meta_tab$Gender)
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
meta_tab$Gender 2 0.1864 0.093197 1.0182 0.0357 0.381
Residuals 55 5.0341 0.091529 0.9643
Total 57 5.2205 1.0000
> adonis(as.dist(d5) ~ meta_tab$milktype)
Call:
adonis(formula = as.dist(d5) ~ meta_tab$milktype)
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
meta_tab$milktype 3 0.4019 0.133976 1.5014 0.07699 0.06 .
Residuals 54 4.8186 0.089233 0.92301
Total 57 5.2205 1.00000
---
Signif. codes: 0 '***’ 0.001 '**’ 0.01 '*’ 0.05 '.’ 0.1 ' ’ 1

第六步:使用PermanovaG()函數(shù)進(jìn)行PERMANOVA對多種UniFrac測量進(jìn)行比較

結(jié)合多個(gè)距離矩陣可以增加假設(shè)檢驗(yàn)的能力,。下面的R代碼使用d(0),,d(0.5),d(1)和函數(shù)PermanovaG()進(jìn)行排列多元方差分析:

> PermanovaG(unifracs[, , c("d_0", "d_0.5", "d_1")] ~ meta_tab$Gestation)
$aov.tab
F.Model p.value
meta_tab$Gestation 1.262764 0.263

9.6  概括和討論

在本章中,我們介紹了多元群落微生物組數(shù)據(jù)的假設(shè)檢驗(yàn)及其在R中的逐步實(shí)現(xiàn),。我們用來說明的數(shù)據(jù)來自我們自己的研究(Jin et al. 2015;王等,?;蚬_提供。讀者可以使用本章提供的R代碼和解釋來分析他們自己的微生物組數(shù)據(jù),。

beta多樣性的多元群落分析有兩種方法:生態(tài)文獻(xiàn)中的“原始數(shù)據(jù)法”與“距離(Mantel)”方法,,PERMANOVA,冗余分析(RDA),,基于距離的RDA(db-RDA),,主坐標(biāo)的規(guī)范化分析屬于基于原始數(shù)據(jù)的方法,而Mantel測試是典型的基于距離的方法,。ANOSIM是一個(gè)基于兩個(gè)距離矩陣之間標(biāo)準(zhǔn)化秩相關(guān)的修正版本的Mantel測試,。因此,它也屬于基于距離的方法,。Mantel試驗(yàn)的應(yīng)用領(lǐng)域也適用于相似度分析(ANOSIM),。

問題是哪一種方法更適合分析多元群落微生物組數(shù)據(jù)?在生態(tài)文學(xué)中,關(guān)于這個(gè)問題的討論很多,。(McArdle and Anderson 2001; Legendre et al. 2005; Legendre2007; Laliberte 2008; Legendre et al. 2008; Pélissier et al. 2008; Tuomisto and Ruokolainen 2008; Anderson et al. 2011; Legendre and Legendre 2012)

基于距離的方法存在以下問題:(1)它們不能正確劃分?jǐn)?shù)據(jù)中的變化,,不能提供正確的第一類錯(cuò)誤率,因此不適合分析beta多樣性;(2) 正如 Anderson 等人認(rèn)為,,更有問題的是,,使用劃分方法對驅(qū)動(dòng)beta多樣性模式的底層過程的相對重要性進(jìn)行直接推斷(Anderson等,2011年),。例如,,ANOSIM和Mantel測試已被證明不適合測試關(guān)于元數(shù)據(jù)變化的假設(shè)檢驗(yàn),因此,,建議在分析beta多樣性時(shí),,應(yīng)限于分析beta多樣性的變化,而不是分析beta多樣性,?;谠紨?shù)據(jù)的替代方法被認(rèn)為為分析beta多樣性提供了更合適和更強(qiáng)大的工具。然而,,基于距離的方法確實(shí)有一些優(yōu)點(diǎn),,例如,Mantel測試被認(rèn)為是一種有效的方法來關(guān)聯(lián)兩個(gè)距離矩陣(Anderson et al.2011)并提供了更大的靈活性,,允許使用其他類型的距離函數(shù),,如Jaccard或Steinhaus/Bray-Curtis (Legendre et al. 2005),它可能適合于關(guān)于群點(diǎn)(樣本)間beta多樣性變化的假設(shè) (Legendre et al. 2005),。與簡單的Mantel試驗(yàn)相比,,生態(tài)學(xué)研究也指出,部分Mantel方法的解釋在一段時(shí)間內(nèi)存在問題.(e.g., Legendre et al. 2005; Anderson et al. 2011).

關(guān)于vegan包中的函數(shù)的使用,這個(gè)包的作者推薦adonis()而不是mrpp()和anosim(),。原因在于函數(shù)adonis()允許用連續(xù)和/或分類預(yù)測器解釋beta多樣性中的方差進(jìn)行類似方差的測試,。然而,mrpp()和anosim()都只處理分類預(yù)測因子,,它們不如adonis()強(qiáng)大,。(Oksanen et al. 2016).

到目前為止,關(guān)于多元群落分析的討論大多局限于生態(tài)學(xué)等相關(guān)研究領(lǐng)域,,而在微生物組領(lǐng)域則沒有,。然而,微生物組學(xué)研究已經(jīng)采用了這些方法和方法,,這些方法的優(yōu)點(diǎn)和局限性也適用于微生物組學(xué)研究,。

根際互作生物學(xué)研究室 簡介


根際互作生物學(xué)研究室是沈其榮教授土壤微生物與有機(jī)肥團(tuán)隊(duì)下的一個(gè)關(guān)注于根際互作的研究小組。本小組由袁軍副教授帶領(lǐng),,主要關(guān)注:1.植物和微生物互作在抗病過程中的作用,;2 環(huán)境微生物大數(shù)據(jù)整合研究;3 環(huán)境代謝組及其與微生物過程研究體系開發(fā)和應(yīng)用,。團(tuán)隊(duì)在過去三年中在 isme J,, Microbiome, PCE,,SBB,,Horticulture Research等期刊上發(fā)表了多篇文章。歡迎關(guān)注 微生信生物 公眾號對本研究小組進(jìn)行了解,。

團(tuán)隊(duì)工作及其成果 (點(diǎn)擊查看)

了解 交流 合作

  • 團(tuán)隊(duì)成員郵箱 袁軍:[email protected],;文濤:[email protected]

  • 團(tuán)隊(duì)公眾號:微生信生物 添加主編微信,或者后臺留言,;

    加主編微信 加入群聊

關(guān)于微生信生物 你想要的都在這里

微生信生物

References Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology 26: 32–46. Anderson, M.J., T.O. Crist, et al. 2011. Navigating the multiple meanings of beta diversity: A roadmap for the practicing ecologist. Ecology Letters 14 (1): 19–28. Benjamini, Y., and Y. Hochberg. 1995. Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B 57: 289–300. Benjamini, Y., and D. Yekutieli. 2001. The control of the false discovery rate in multiple testing under dependency. Annals of Statistics 29: 1165–1188. Bonferroni, C.E. 1936. Teoria statistica delle classi e calcolo delle probabilitàby. Pubblicazioni del R Istituto Superiore di Scienze Economiche e Commerciali di Firenze 8: 3–62 Key: citeulike:1778138. Carcer, D.A., S.E. Denman, et al. 2011. Evaluation of subsampling-based normalization strategies for tagged high-throughput sequencing data sets from gut microbiomes. Applied and Environment Microbiology 77 (24): 8795–8798. Chang, Q., Y. Luan, et al. 2011. Variance adjusted weighted UniFrac: A powerful beta diversity measure for comparing communities based on phylogeny. BMC Bioinformatics 12 (1): 118. Chen, J. 2012. GUniFrac: Generalized UniFrac distances. R package version 1.0. https://CRAN.Rproject.org/package=GUniFrac. Chen, J., K. Bittinger, et al. 2012. Associating microbiome composition with environmental covariates using generalized UniFrac distances. Bioinformatics 28 (16): 2106–2113. Clarke, K.R. 1993. Non-parametric multivariate analysis of changes in community structure. Australian Journal of Ecology 18: 117–143. Hochberg, Y. 1988. A sharper Bonferroni procedure for multiple tests of significance. Biometrika 75: 800–803. Holm, S. 1979. A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics 6: 65–70. Hommel, G. 1988. A stagewise rejective multiple test procedure based on a modified Bonferroni test. Biometrika 75: 383–386. Jin, D., S. Wu, et al. 2015. Lack of Vitamin D receptor causes dysbiosis and changes the functions of the murine intestinal microbiome. Clinical Therapeutics 37 (5): 996–1009, e1007. Laliberte, E. 2008. Analyzing or explaining beta diversity? comment. Ecology 89 (11): 3232– 3237. Legendre, P. 2007. Studying beta diversity: Ecological variation partitioning by multiple regression and canonical analysis. Journal of Plant Ecology 1 (1): 3–8. Legendre, P., and L. Legendre. 2012. Numerical ecology. Amsterdam: Elsevier Science BV. Legendre, P., and M.J. Fortin. 1989. Spatial pattern and ecological analysis. Vegetatio 80 (2): 107–138. Legendre, P., D. Borcard, et al. 2005. Analyzing beta diversity: Partitioning the spatial variation of community composition data. Ecological Monographs 75 (4): 435–450. Legendre, P., D. Borcard, et al. 2008. Analyzing or explaining beta diversity? comment. Ecology 89 (11): 3238–3244. 9.6 Summary and Discussion 329 Linnenbrink, M., J. Wang, et al. 2013. The role of biogeography in shaping diversity of the intestinal microbiota in house mice. Molecular Ecology 22 (7): 1904–1916. Lozupone, C., and R. Knight. 2005. UniFrac: A new phylogenetic method for comparing microbial communities. Applied and Environmental Microbiology 71 (12): 8228–8235. Lozupone, C.A., M. Hamady, et al. 2007. Quantitative and qualitative beta diversity measures lead to different insights into factors that structure microbial communities. Applied and Environment Microbiology 73 (5): 1576–1585. Mantel, N. 1967. The detection of disease clustering and a generalized regression approach. Cancer Research 27: 209–220. Mantel, N., and R.S. Valand. 1970. A technique of nonparametric multivariate analysis. Biometrics 26 (3): 547–558. McArdle, B.H., and M.J. Anderson. 2001. Fitting multivariate models to community data: A comment on distance-based redundancy analysis. Ecology 82 (1): 290–297. Mielke, P.W. 1984. Meteorological applications of permutation techniques based on distance functions. In Handbook of statistics, vol. 4, ed. P.R. Krishnaiah and P.K. Sen, 813–830. Amsterdam, North-Holland: Elsevier Science Publishers. Mielke Jr., P.W. 1991. The application of multivariate permutation methods based on distance functions in the earth sciences. Earth-Science Reviews 31: 55–71. Oksanen, J., F. Guillaume Blanchet, et al. 2016. Vegan: Community ecology package. R package version 2.4-1. http://CRAN.R-project.org/package=vegan. Pélissier, R., P. Couteron, et al. 2008. Analyzing or explaining beta diversity? comment. Ecology 89 (11): 3227–3232. Smouse, P.E., J.C. Long, et al. 1986. Multiple regression and correlation extensions of the mantel test of matrix correspondence. Systematic Zoology 35 (4): 627–632. Tuomisto, H., and K. Ruokolainen. 2006. Analyzing or explaining beta diversity? Understanding the targets of different methods of analysis. Ecology 87: 2697–2708. Tuomisto, H., and K. Ruokolainen. 2008. Analyzing or explaining beta diversity: Reply. Ecology 89: 3244–3256. Urbaniak, C., M. Angelini, et al. 2016. Human milk microbiota profiles in relation to birthing method, gestation and infant gender. Microbiome 4 (1): 1.

    轉(zhuǎn)藏 分享 獻(xiàn)花(0

    0條評論

    發(fā)表

    請遵守用戶 評論公約

    類似文章