微生物組數據分析方法（二）

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

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

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

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

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

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

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

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

常用單因素方差分析（one-way anova）來(lái)進(jìn)行比較，其適用于只有一個(gè)因素（自變量）的處理。在觀(guān)測變量總離差平方和中，如果組間離差平方和所占比例較大，則說(shuō)明觀(guān)測變量的變動(dòng)主要是由控制變量引起的，可以主要由控制變量來(lái)解釋?zhuān)刂谱兞拷o觀(guān)測變量帶來(lái)了顯著(zhù)影響；反之，如果組間離差平方和所占比例小，則說(shuō)明觀(guān)測變量的變動(dòng)不是主要由控制變量引起的，不可以主要由控制變量來(lái)解釋?zhuān)刂谱兞康牟煌經(jīng)]有給觀(guān)測變量帶來(lái)顯著(zhù)影響，觀(guān)測變量值的變動(dòng)是由隨機變量因素引起的。

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

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

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

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

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

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

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

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

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

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

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

微生物組數據分析方法之：功能基因預測

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

微生物組數據分析方法之：microPITA分析

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

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

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

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

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

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

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

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

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

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

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

微生物組數據分析方法之：溯源分析

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

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

參考文獻：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)系我們，我們將免費為您設計文章思路方案