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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

BugBase：是一種預測復雜微生物組內功能途徑的生物水平覆蓋以及生物可解釋表型的方法。BugBase首先通過預測的16S拷貝數對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)學中的概念，其涉及到一類物種（不管親緣關系相近與否），這些物種能通過相似的途徑利用同類環(huán)境資源。

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

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

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

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

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

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

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

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

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

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

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

微生物組數據分析方法之：相關性與關聯分析

RDA/CCA分析：RDA(Redundancy analysis)/CCA(Canonical Correspondence analysis)是基于對應分析發(fā)展的一種排序方法，RDA分析基于線性模型，CCA分析基于單峰模型，主要用來反映菌群或樣品與環(huá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距離等，而不再僅僅局限于歐氏距離）進行RDA排序，因此db-RDA在生態(tài)學統計分析中被廣泛使用。當然，若是我們直接計算樣方群落間的歐式距離矩陣并將其輸入至db-RDA中計算時，也將會得到和常規(guī)的RDA一致的結果。

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

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

網絡通常由邊和節(jié)點構成，一條邊由兩個節(jié)點連接而成。網絡中節(jié)點具有很多重要的性質，包括度、聚類系數、緊密中心性、中介中心性、Zi（within-module connectivity）與 Pi（among-module connectivity）等。通常度、聚類系數、緊密中心性、中介中心性值越大，說明節(jié)點重要性越高；此外通過Zi和Pi值可以將節(jié)點分為四類：Peripheral nodes（zi?≤?2.5, Pi?≤?0.62，僅有少量的邊并且通常只與模塊內部的節(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é)點高度連接又能連接不同的模塊)，通過這些性質可以說明網絡中節(jié)點的重要性。

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

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

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

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

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

目標樣本為Sink，微生物污染源或來源的樣品為Source；基于貝葉斯算法，探究目標樣本（Sink）中微生物污染源或來源（Source）的分析。根據Source樣本和Sink樣本的群落結構分布，來預測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.

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