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

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

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

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

(2) 隨機誤差，如測量誤差造成的差異或個體間的差異，稱為組內(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。另一種情況是處理確實有作用，組間均方是由于誤差與不同處理共同導致的結(jié)果，即各樣本來自不同總體，那么，MSb>>MSw(遠遠大于)。MSb/MSw比值構(gòu)成F分布，用F值與其臨界值比較，推斷各樣本是否來自相同的總體。

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

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

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

該分析針對不同樣本會采用不同分析方法：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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