微生物組數(shù)據(jù)分析方法之：組間差異顯著性分析

上篇介紹了關(guān)于群落結(jié)構(gòu)整體間的差異分析，那么針對不同的處理，不同的環(huán)境條件下，究竟是哪些微生物對環(huán)境變化更加敏感，在處理組間表現(xiàn)出更顯著的差異，而這些微生物又是如何響應(yīng)環(huán)境變化的，是否可以作為指示物種用于后續(xù)的實(shí)驗(yàn)研究當(dāng)中，這些就涉及到我們今天要講的組間差異顯著性分析，一起來看下：

ANOVA分析：方差分析(Analysis of Variance，簡稱ANOVA)，又稱“變異數(shù)分析”或“F檢驗(yàn)”， 用于兩個(gè)及兩個(gè)以上樣本均數(shù)差別的顯著性檢驗(yàn)。它的基本思想是將測量數(shù)據(jù)的總變異（即總方差）按照變異來源分為處理（組間）效應(yīng)和誤差（組內(nèi)）效應(yīng)，通過分析研究不同來源的變異對總變異的貢獻(xiàn)大小，從而確定可控因素對研究結(jié)果影響力的大小。

(1) 實(shí)驗(yàn)條件，即不同的處理造成的差異，稱為組間差異，用變量在各組的均值與總均值之偏差平方和的總和表示，記作SSb，組間自由度dfb。

(2) 隨機(jī)誤差，如測量誤差造成的差異或個(gè)體間的差異，稱為組內(nèi)差異，用變量在各組的均值與該組內(nèi)變量值之偏差平方和的總和表示，記作SSw，組內(nèi)自由度dfw。

總偏差平方和 SSt = SSb + SSw。

組內(nèi)SSw、組間SSb除以各自的自由度(組內(nèi)dfw =n-m，組間dfb=m-1，其中n為樣本總數(shù)，m為組數(shù))，得到其均方MSw和MSb。

一種情況是處理沒有作用，即各組樣本均來自同一總體，MSb/MSw≈1。另一種情況是處理確實(shí)有作用，組間均方是由于誤差與不同處理共同導(dǎo)致的結(jié)果，即各樣本來自不同總體，那么，MSb>>MSw(遠(yuǎn)遠(yuǎn)大于)。MSb/MSw比值構(gòu)成F分布，用F值與其臨界值比較，推斷各樣本是否來自相同的總體。

常用單因素方差分析（one-way anova）來進(jìn)行比較，其適用于只有一個(gè)因素（自變量）的處理。在觀測變量總離差平方和中，如果組間離差平方和所占比例較大，則說明觀測變量的變動(dòng)主要是由控制變量引起的，可以主要由控制變量來解釋，控制變量給觀測變量帶來了顯著影響；反之，如果組間離差平方和所占比例小，則說明觀測變量的變動(dòng)不是主要由控制變量引起的，不可以主要由控制變量來解釋，控制變量的不同水平?jīng)]有給觀測變量帶來顯著影響，觀測變量值的變動(dòng)是由隨機(jī)變量因素引起的。

秩和檢驗(yàn)：秩和檢驗(yàn)是非參數(shù)檢驗(yàn)的一種方法。其中“秩”又稱等級(jí)、即上述次序號(hào)的和稱“秩和”，秩和檢驗(yàn)就是用秩和作為統(tǒng)計(jì)量進(jìn)行假設(shè)檢驗(yàn)的方法。

當(dāng)在實(shí)踐中遇到總體分布類型未知；或者總體分布類型已知但不符合正態(tài)分布；或者某些變量可能無法jing q測量；以及方差不齊時(shí)，即可采用秩和檢驗(yàn)進(jìn)行分析。

該分析針對不同樣本會(huì)采用不同分析方法：

1）兩組獨(dú)立樣本 — 小樣本采用Wilcoxon秩和檢驗(yàn) （Wilcoxon rank-sum test），大樣本采用曼-惠特尼 U 檢驗(yàn)（Mann–Whitney U test），樣本數(shù)小于等于20為小樣本，樣本數(shù)大于20為大樣本；

2）多組獨(dú)立樣本 — Kruskal-Wallis秩和檢驗(yàn)，兩兩分析可采用Nemenyi法檢驗(yàn)，Kruskal-Wallis秩和檢驗(yàn)要求每組樣本數(shù)不得小于5個(gè)。

該分析可以對兩組或多組樣品的物種、基因或者功能進(jìn)行顯著性差異分析，并對 p值進(jìn)行FDR校正。

LEfSe分析：LEfSe (Line Discriminant Analysis (LDA) Effect Size)分析是一種將非參數(shù)的Kruskal-Wallis以及Wilcoxon秩和檢驗(yàn)，與線性判別分析(Linear discriminant analysis，LDA)效應(yīng)量(Effect size)相結(jié)合的分析手段, 能夠在不同組間尋找具有統(tǒng)計(jì)學(xué)差異的Biomarker，其要求組內(nèi)樣本數(shù)≥3。

首先使用non-parametric factorial Kruskal-Wallis (KW) sum-rank test（非參數(shù)因子克魯斯卡爾—沃利斯和秩驗(yàn)檢）檢測具有顯著豐度差異特征，并找到與豐度有顯著性差異的類群，然后采用線性判別分析（LDA）來估算每個(gè)組分（物種）豐度對差異效果影響的大小。

Metastats分析：用于尋找組間差異物種，從算法原理上來說，Metastats實(shí)際上是非參數(shù)多重檢驗(yàn)和p值校正的整合，由于Metastats采用的非參數(shù)t檢驗(yàn)，因此只能分析兩個(gè)分組間差異。對組間的物種豐度數(shù)據(jù)進(jìn)行 T 檢驗(yàn)得到 p 值，并對 p 值進(jìn)行校正得到 q 值；最后根據(jù) p 值（或 q 值）篩選出導(dǎo)致兩組樣品組成差異的物種，進(jìn)而進(jìn)行可視化繪圖，得到相應(yīng)的物種差異柱狀圖，或者箱線圖展示。

Metagenomeseq分析：Metagenomseq方法是近幾年出現(xiàn)的新的兩組間檢驗(yàn)的方法，用normalization實(shí)現(xiàn)分類注釋時(shí)的biases處理，同時(shí)用零膨脹高斯分布（zero-flated Gaussian distribution）處理了測序深度所帶來的影響，在此基礎(chǔ)上，利用線性模型找到存在的差異所在。用于研究兩組樣品間微生物群落豐度的差異。針對其分析結(jié)果，也可以通過可視化繪圖，得到相應(yīng)的物種差異柱狀圖，或者箱線圖等。

隨機(jī)森林分析：Random forest 分析，隨機(jī)森林指的是利用多棵樹對樣本進(jìn)行訓(xùn)練并預(yù)測的一種分類器，是機(jī)器學(xué)習(xí)的一種。利用隨機(jī)森林分析可以幫我們找到數(shù)據(jù)中重要的特征：在微生物組中就是哪些微生物組成可以作為特征進(jìn)行各種預(yù)測；另外，隨機(jī)森林可以利用這些重要特征構(gòu)建模型用于分類預(yù)測（分類變量）或者回歸預(yù)測（數(shù)值型變量）；隨機(jī)森林分析建議分組內(nèi)樣品至少30個(gè)以上，樣品數(shù)太少影響準(zhǔn)確性。 隨機(jī)森林分析，一般將數(shù)據(jù)劃分為訓(xùn)練數(shù)據(jù)集（train data）和測試數(shù)據(jù)集（test or validate data），訓(xùn)練數(shù)據(jù)集用于構(gòu)建預(yù)測模型，測試集用于查看模型的準(zhǔn)確性。

微生物組數(shù)據(jù)分析方法之：功能基因預(yù)測

PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States)是一款基于樣本中的標(biāo)記基因序列豐度來預(yù)測樣本功能豐度的軟件。這里的功能是指基因家族，如KEGG同源基因、EC酶分類號(hào)等，其原理如下圖：

PICRUSt 2 將特征序列（16S rRNA）與微生物基因組數(shù)據(jù)庫（IMG，Integrated Microbial Genomes）參考序列對齊（aliign）構(gòu)建進(jìn)化樹，找到特征序列的“最近物種”，并根據(jù)已知物種的基因種類和豐度信息預(yù)測未知物種的基因信息，從而結(jié)合基因的KEGG通路信息預(yù)測整個(gè)群落的通路情況。通過組間差異分析，得到在組間具有顯著差異的功能，最后利用Stamp對結(jié)果進(jìn)行差異可視化。

TAX4FUN2： Tax4Fun2可基于16S rRNA基因序列快速預(yù)測原核生物的功能譜和功能冗余。通過合并用戶定義的、特定于棲息地的基因組信息，可以顯著提高預(yù)測功能圖譜的準(zhǔn)確性。不再局限于僅SILVA的特定版本注釋的OTU豐度表，允許直接以O(shè)TU代表序列作為輸入，通過與指定參考數(shù)據(jù)庫的比對實(shí)現(xiàn)物種注釋。除了Tax4Fun2提供的已構(gòu)建好的參考集（相比之前大幅擴(kuò)大），也允許我們提供自定義的參考集，使用非常靈活。

FAPROTAX：較適用于對環(huán)境樣本（如海洋、湖泊等）的生物地球化學(xué)循環(huán)過程（特別是碳、氫、氮、磷、硫等元素循環(huán)）等生態(tài)功能進(jìn)行預(yù)測。FAPROTAX 是根據(jù)已發(fā)表的文獻(xiàn)手動(dòng)構(gòu)建的數(shù)據(jù)庫，它把原核微生物的分類和代謝等功能對應(yīng)起來。因其基于已發(fā)表驗(yàn)證的可培養(yǎng)菌文獻(xiàn)，其預(yù)測準(zhǔn)確度較好，但相比于上述PICRUSt2和Tax4Fun來說預(yù)測的覆蓋度可能會(huì)降低。

BugBase：是一種預(yù)測復(fù)雜微生物組內(nèi)功能途徑的生物水平覆蓋以及生物可解釋表型的方法。BugBase首先通過預(yù)測的16S拷貝數(shù)對OTU進(jìn)行歸一化，然后使用提供的預(yù)先計(jì)算的文件預(yù)測微生物表型。包括以下九個(gè)方面：Aerobic 好氧、Anaerobic 厭氧、Contains Mobile Elements 移動(dòng)元件、Facultatively Anaerobic 兼性厭氧、Forms Biofilms ?生物膜形成

Gram Negative 革蘭氏陰性、Gram Positive 革蘭氏陽性、Potentially Pathogenic 潛在致病性、Stress Tolerant 脅迫耐受。

FUNGuild：FUNGuild（ Fungi Functional Guild）是一款通過微生態(tài) guild 對真菌群落進(jìn)行分類分析的工具，guild 是微生態(tài)學(xué)中的概念，其涉及到一類物種（不管親緣關(guān)系相近與否），這些物種能通過相似的途徑利用同類環(huán)境資源。

FUNGuild 基于目前已發(fā)表的文獻(xiàn)或權(quán)威網(wǎng)站數(shù)據(jù)，首先根據(jù)營養(yǎng)方式（trophicMode分類系統(tǒng)）將真菌分為三大類：

病理營養(yǎng)型（pathotroph）：通過損害宿主細(xì)胞而獲取營養(yǎng)（包括吞噬型真菌 phagotrophs）；

共生營養(yǎng)型（symbiotroph）：通過與宿主細(xì)胞交換資源來獲取營養(yǎng)；

腐生營養(yǎng)型（saprotroph）：通過降解死亡的宿主細(xì)胞來獲取營養(yǎng)。

基于 3 大營養(yǎng)方式，又進(jìn)一步細(xì)分為若干個(gè) guilds，用來對trophicMode 分類系統(tǒng)的補(bǔ)充：

在Pathotroph 下，細(xì)分為 Animal Pathogen ： 動(dòng)物病原菌、Plant Pathogen：植物病原菌、Fungal Parasite：真菌寄生菌、Lichen Parasite：地衣寄生菌、Bryophyte Parasite：苔蘚植物寄生菌、Clavicipitaceous Endophyte：內(nèi)生真菌。

在Symbiotroph 下，細(xì)分為Ectomycorrhizal：外生菌根、Ericoid Mycorrhizal：杜鵑花類菌根、Lichenized：地衣共生菌等。

在Saprotroph 下，細(xì)分成Dung Saprotroph：排泄物腐生菌（如糞便）、Leaf Saprotroph：葉子腐生菌、Plant Saprotroph：植物腐生菌 （生長環(huán)境多腐敗的植物）、Soil Saprotroph：土壤腐生菌、Wood Saprotroph：木質(zhì)腐生菌。

基于我們測序得到的真菌序列，就可以進(jìn)行Guild 的功能注釋。

微生物組數(shù)據(jù)分析方法之：microPITA分析

MicroPITA：在做完16S、18S或ITS等微生物多樣性研究后，我們常常還會(huì)想進(jìn)一步了解微生物群落的功能。通常情況下，會(huì)采用宏基因組、宏轉(zhuǎn)錄組或宏代謝組等方法深入分析，但相對于擴(kuò)增子測序，宏基因組等測序手段的價(jià)格還是相對較高，因此需要從已測完的樣本中再挑選合適的樣本進(jìn)行宏基因組測序?？衫胢icroPITA進(jìn)行樣品預(yù)測，挑選出合適的樣品。該分析是基于大量微生物多樣性的數(shù)據(jù)，根據(jù)不同指標(biāo)篩選出代表性樣本，以便于開展開展后續(xù)研究。

微生物組數(shù)據(jù)分析方法之：相關(guān)性與關(guān)聯(lián)分析

RDA/CCA分析：RDA(Redundancy analysis)/CCA(Canonical Correspondence analysis)是基于對應(yīng)分析發(fā)展的一種排序方法，RDA分析基于線性模型，CCA分析基于單峰模型，主要用來反映菌群或樣品與環(huán)境因子之間的關(guān)系。根據(jù)物種分布變化，選擇最*的分析模型。

RDA或CCA模型的選擇原則：先用species-sample豐度數(shù)據(jù)做DCA分析，看分析結(jié)果中Lengths of gradient 的第一軸的大小，如果大于4.0，就應(yīng)該選CCA，如果3.0-4.0之間，選RDA和CCA均可，如果小于3.0，RDA的結(jié)果要好于CCA。

db-RDA分析：基于距離的冗余分析（db-RDA），db-RDA將主坐標(biāo)分析（PCoA）計(jì)算的樣方得分矩陣應(yīng)用在RDA中，其好處是可以基于任意一種距離測度（例如Bray-curtis距離等，而不再僅僅局限于歐氏距離）進(jìn)行RDA排序，因此db-RDA在生態(tài)學(xué)統(tǒng)計(jì)分析中被廣泛使用。當(dāng)然，若是我們直接計(jì)算樣方群落間的歐式距離矩陣并將其輸入至db-RDA中計(jì)算時(shí)，也將會(huì)得到和常規(guī)的RDA一致的結(jié)果。

Mantel Test分析：Mantel test 是檢驗(yàn)兩個(gè)矩陣相關(guān)關(guān)系的非參數(shù)統(tǒng)計(jì)方法，多用在生態(tài)學(xué)上，用來檢驗(yàn)群落距離矩陣(如 Bray-Curtis distance matrix)和環(huán)境變量距離矩陣（即環(huán)境因子指標(biāo)）之間的相關(guān)性。其可以得到特征與環(huán)境因子之間的mantel test的r值和p值。

相關(guān)性網(wǎng)絡(luò)分析：網(wǎng)絡(luò)圖是相關(guān)性分析的一種表現(xiàn)形式，根據(jù)各個(gè)物種在各個(gè)樣品中的豐度以及變化情況，進(jìn)行斯皮爾曼(Spearman )秩相關(guān)分析并篩選相應(yīng)條件下的數(shù)據(jù)構(gòu)建相關(guān)性網(wǎng)絡(luò)?；诰W(wǎng)絡(luò)圖的分析，可以獲得物種在環(huán)境樣本中的共存關(guān)系，得到物種在同一環(huán)境下的相互作用的情況及重要的模式信息，進(jìn)一步解釋樣本間表型差異的形成機(jī)制。

網(wǎng)絡(luò)通常由邊和節(jié)點(diǎn)構(gòu)成，一條邊由兩個(gè)節(jié)點(diǎn)連接而成。網(wǎng)絡(luò)中節(jié)點(diǎn)具有很多重要的性質(zhì)，包括度、聚類系數(shù)、緊密中心性、中介中心性、Zi（within-module connectivity）與 Pi（among-module connectivity）等。通常度、聚類系數(shù)、緊密中心性、中介中心性值越大，說明節(jié)點(diǎn)重要性越高；此外通過Zi和Pi值可以將節(jié)點(diǎn)分為四類：Peripheral nodes（zi?≤?2.5, Pi?≤?0.62，僅有少量的邊并且通常只與模塊內(nèi)部的節(jié)點(diǎn)相連）、Connectors (zi?≤?2.5, Pi?>?0.62，通常連接不同的模塊)、Module hubs (zi?>?2.5, Pi?≤?0.62，與自身所在模塊中的許多節(jié)點(diǎn)高度連接)和Network hubs (zi?>?2.5, Pi?>?0.62，即能與自身所在模塊中的許多節(jié)點(diǎn)高度連接又能連接不同的模塊)，通過這些性質(zhì)可以說明網(wǎng)絡(luò)中節(jié)點(diǎn)的重要性。

網(wǎng)絡(luò)自身也具有很多不同的特性，例如節(jié)點(diǎn)數(shù)據(jù)量、邊數(shù)量、模塊性與模塊數(shù)量、網(wǎng)絡(luò)直徑與密度、平均路徑與平均聚類系數(shù)等，共同體現(xiàn)網(wǎng)絡(luò)屬性。

VIF(方差膨脹因子)分析：影響菌群的環(huán)境因子因素很多，但是有些環(huán)境因子間存在較強(qiáng)的多重共線性相關(guān)關(guān)系，會(huì)影響后續(xù)分析，因此在環(huán)境因子關(guān)聯(lián)分析前，可對環(huán)境因子進(jìn)行篩選，保留共線性較小的環(huán)境因子用于后續(xù)分析。方差膨脹因子（Variance Inflation Factor，VIF）分析時(shí)目前常用的環(huán)境因子篩選方法，VIF是衡量多元線性回歸模型中復(fù)雜(多重)共線性嚴(yán)重程度的一種度量。它表示回歸系數(shù)估計(jì)量的方差與假設(shè)自變量間不線性相關(guān)時(shí)方差相比的比值。其計(jì)算公式為：VIFi =1/（1-Ri^2），其中Ri^2代表模型中與其它自變量相關(guān)的第i個(gè)自變量的方差比例，用于衡量第i個(gè)自變量與其它自變量間的共線性關(guān)系。VIF越大，顯示共線性越嚴(yán)重。經(jīng)驗(yàn)判斷方法表明：當(dāng)VIF大于0且小于10，不存在多重共線性；當(dāng)VIF大約等于10且小于100，存在較強(qiáng)的多重共線性；當(dāng)VIF大于等于100，存在嚴(yán)重多重共線性，因此通常認(rèn)為小于10的環(huán)境因子是無用的環(huán)境因子。

VPA分析(方差分解分析)：VPA，全稱Variance Partitioning Analysis，中文稱為方差分解分析，該分析的目的是確定指定的環(huán)境因子對群落結(jié)構(gòu)變化的解釋比例。通常在CCA/RDA的排序分析方法可以得到所有參與分析的環(huán)境因子對群落變化的解釋比例。使用R語言vegan包(version:2.3-0)，varpart函數(shù)進(jìn)行分析；分析前首先需要對環(huán)境因子進(jìn)行一個(gè)分類，一般支持2-4類，默認(rèn)使用大于等于2個(gè)環(huán)境因子進(jìn)行此分析。

微生物組數(shù)據(jù)分析方法之：溯源分析

SourceTracker 分析：SourceTracker—微生物來源分析，用途是可以識(shí)別相關(guān)各組間來源的分析，如嬰兒的腸道菌群有哪些繼承了母親的腸道菌群、哪些來自陰道菌群、哪些來自皮膚；河流污染物的來源分析、周圍工廠、農(nóng)田、養(yǎng)殖廠對河流污染的貢獻(xiàn)和來源追溯等?？梢怨烙?jì)末知群體微生物來源和組成，有用的是可根據(jù)環(huán)境樣品來分類微生物的來源。

目標(biāo)樣本為Sink，微生物污染源或來源的樣品為Source；基于貝葉斯算法，探究目標(biāo)樣本（Sink）中微生物污染源或來源（Source）的分析。根據(jù)Source樣本和Sink樣本的群落結(jié)構(gòu)分布，來預(yù)測Sink樣本中來源于各Source樣本的組成比例。

參考文獻(xiàn)：Knights D, Kuczynski J, Charlson ES, et al. Bayesian community-wide culture-independent microbial source tracking [J]. Nature Methods, 2011, 8(9): 761.

如果您對以上分析方案感興趣，歡迎點(diǎn)擊下方按鈕聯(lián)系我們，我們將免費(fèi)為您設(shè)計(jì)文章思路方案