樊 哲 呂嘉麗 張 帥 谷冰冰 張 濤△

【提 要】 目的 微生物組研究中的典型數(shù)據(jù)是由操作分類單元(OTU)計數(shù)組成,這些數(shù)據(jù)具有零膨脹、過離散、偏態(tài)等特點,常常被研究者忽視。本文旨在介紹微生物組學(xué)數(shù)據(jù)關(guān)聯(lián)分析方法,并提出相關(guān)研究策略。方法 根據(jù)微生物組學(xué)數(shù)據(jù)特點,本文擬從線性相關(guān)分析和非線性相關(guān)分析的角度,對近年來國內(nèi)外研究者提出的微生物組學(xué)數(shù)據(jù)關(guān)聯(lián)分析的統(tǒng)計分析方法展開討論。結(jié)果 根據(jù)各方法的核心思想及優(yōu)缺點,總結(jié)統(tǒng)計分析策略。結(jié)論 運用本文提出的統(tǒng)計分析策略,可有效識別微生物同其他研究指標(biāo)間多種函數(shù)形式的關(guān)聯(lián)。為探討微生物對人體健康影響,研究微生物的作用機制提供線索。

近年來,隨著高通量測序技術(shù)及統(tǒng)計分析方法的發(fā)展,微生物組學(xué)聯(lián)合代謝組、蛋白組及其他組學(xué)的多組學(xué)關(guān)聯(lián)研究已成為系統(tǒng)生物學(xué)研究新趨勢[1]。從系統(tǒng)生物學(xué)角度來看,整合多個層面組學(xué)數(shù)據(jù),構(gòu)建多組學(xué)關(guān)聯(lián)網(wǎng)絡(luò),能更充分理解各分子間的調(diào)控及因果關(guān)系,為探索疾病深層發(fā)生發(fā)展機制提供新思路。組學(xué)數(shù)據(jù)通常具有高維的特點,在探索分析的階段,可以通過關(guān)聯(lián)分析的統(tǒng)計方法,過濾出有統(tǒng)計學(xué)意義的變量,從而為后續(xù)的機制分析篩選出更小范圍的目標(biāo)變量。然而,目前多數(shù)研究在微生物組關(guān)聯(lián)分析中往往忽視了微生物組學(xué)數(shù)據(jù)特點及各相關(guān)性分析方法的適用條件,導(dǎo)致研究結(jié)果產(chǎn)生偏差,出現(xiàn)關(guān)聯(lián)結(jié)論與生物學(xué)結(jié)論不一致的問題[2-5]。本文擬對近年來國內(nèi)外研究者提出的微生物組學(xué)關(guān)聯(lián)分析方法進行介紹,并系統(tǒng)地總結(jié)各個方法的核心思想及優(yōu)缺點,提出微生物組學(xué)數(shù)據(jù)關(guān)聯(lián)分析的統(tǒng)計分析策略。

微生物組學(xué)數(shù)據(jù)特點

微生物組學(xué)數(shù)據(jù)通常來自16S測序或宏基因組測序。16S測序是在提取微生物DNA后,對微生物16S rRNA基因高變異區(qū)域進行聚合酶鏈?zhǔn)椒磻?yīng) (polymerase chain reaction,PCR) 擴增和測序。而宏基因組測序是對研究樣本中全部微生物的總DNA進行高通量測序,并且宏基因組測序有更深的測序深度,能鑒定到種水平甚至菌株水平的微生物[6]。處理后的序列在一定的相似度水平上聚類為操作分類單元 (operational taxonomic units,OTU),通常將相似水平大于97%的OTU納入后續(xù)生物信息學(xué)和統(tǒng)計學(xué)分析。

真實的微生物組學(xué)數(shù)據(jù)(如表1所示)具有以下幾個特征:① 稀疏性:觀測單位中,每一個OTU有部分觀測值為零,也被稱為零膨脹現(xiàn)象。② 非線性:微生物豐度與其他組學(xué)數(shù)據(jù)或臨床指標(biāo)存在多種多樣的非線性關(guān)系,而不只是單一的線性相關(guān)。③ 高維特點:微生物組學(xué)數(shù)據(jù)還具有與其他組學(xué)數(shù)據(jù)相同的高維特點,即變量數(shù)大于或遠大于樣本數(shù)的情況。④ 過離散:OTU數(shù)據(jù)存在過離散的現(xiàn)象,即變量的方差大于均值。迄今為止,仍沒有一種方法或策略能夠應(yīng)對微生物數(shù)據(jù)的眾多挑戰(zhàn)。

表1 微生物組學(xué)數(shù)據(jù)分布

微生物組學(xué)關(guān)聯(lián)分析常用統(tǒng)計分析方法

微生物組學(xué)關(guān)聯(lián)分析的研究目的通常包括:①從高維特征中過濾出少量特征用于后續(xù)機制研究。②構(gòu)建多組學(xué)關(guān)聯(lián)網(wǎng)絡(luò)。傳統(tǒng)的關(guān)聯(lián)分析方法僅能發(fā)現(xiàn)兩變量之間的單調(diào)關(guān)系,而微生物同其他水平分子之間可能存在多種函數(shù)形式的非線性關(guān)系。因此,根據(jù)研究目的和相關(guān)關(guān)系的形式,微生物組學(xué)相關(guān)性分析方法可分為以下幾類:①傳統(tǒng)關(guān)聯(lián)分析方法:包括傳統(tǒng)單變量相關(guān)統(tǒng)計方法及其用于微生物組數(shù)據(jù)的擴展方法,僅能檢驗出線性關(guān)系。②零膨脹模型:檢驗變量間的線性關(guān)系,更好地解決了微生物組數(shù)據(jù)偏態(tài),零膨脹和過離散等問題。③基于互信息的方法:檢驗變量間的非線性關(guān)聯(lián)關(guān)系,且無分布要求。④知識驅(qū)動的方法:構(gòu)建關(guān)聯(lián)網(wǎng)絡(luò)。本研究對十種關(guān)聯(lián)分析方法進行簡要的介紹和比較(如表2所示)。

表2 關(guān)聯(lián)分析方法原理及優(yōu)點

1.傳統(tǒng)關(guān)聯(lián)分析方法及其擴展方法

線性相關(guān)分析因具有可解釋性強的優(yōu)勢而經(jīng)常被使用。Pearson相關(guān)以兩變量與各自平均的離差反映兩連續(xù)變量的線性相關(guān)程度[7]。該方法可以定量描述相關(guān)強度與方向,但當(dāng)不滿足正態(tài)性假設(shè)或樣本量小于30時,難以得到可信結(jié)果。Mao等在山羊瘤胃微生物組的多組學(xué)研究中,應(yīng)用單變量相關(guān)方法建立微生物的屬與代謝物之間的Pearson相關(guān)矩陣,發(fā)現(xiàn)瘤胃微生物群落結(jié)構(gòu)的改變與代謝物之間存在明顯的相關(guān)性[8]。然而Pearson方法對變量分布要求嚴(yán)格,且微生物組數(shù)據(jù)通常呈現(xiàn)偏態(tài)或極偏態(tài)分布,不滿足Pearson方法的使用條件,繼續(xù)使用該方法可能導(dǎo)致檢驗效能下降,甚至出現(xiàn)錯誤的結(jié)果。

Friedman等在Pearson相關(guān)分析的基礎(chǔ)上進行了改進,提出了成分?jǐn)?shù)據(jù)稀疏關(guān)聯(lián)(sparse correlations for compositional data,SparCC) 方法用于分析微生物成分?jǐn)?shù)據(jù),根據(jù)觀測值的Dirichlet分布對關(guān)聯(lián)程度進行估計[9]。微生物組中的成分?jǐn)?shù)據(jù)即相對豐度數(shù)據(jù),描述單個物種占整個樣本微生物群落的百分比,某一樣本內(nèi)所有物種的相對豐度相加等于百分之百這一特點可能會導(dǎo)致傳統(tǒng)的關(guān)聯(lián)分析方法得到虛假的關(guān)聯(lián)?；趌asso的成分?jǐn)?shù)據(jù)關(guān)聯(lián)方法(correlation inference for compositional data through lasso,CCLasso)也適用于微生物數(shù)據(jù)關(guān)聯(lián)分析,其原理是在對成分?jǐn)?shù)據(jù)進行對數(shù)轉(zhuǎn)換后,使用基于懲罰函數(shù)的最小二乘法估計相關(guān)系數(shù)[10]。SparCC方法用于成分?jǐn)?shù)據(jù)分析,解決了一般方法結(jié)果中正相關(guān)性被抑制的問題。You等人通過模擬研究發(fā)現(xiàn),相較于Pearson方法,在微生物和代謝物關(guān)聯(lián)分析時使用SparCC和CCLasso兩種方法能更好地控制假陽性錯誤率[11]。這兩種方法作為傳統(tǒng)方法在微生物組學(xué)中的擴展,雖然解決了微生物組數(shù)據(jù)偏態(tài)的問題,但忽視了微生物組中數(shù)據(jù)的零膨脹特點,觀測值中過多的零值可能導(dǎo)致這些方法性能下降。

Spearman秩相關(guān)分析適用于估計連續(xù)或有序變量間的單調(diào)關(guān)聯(lián)研究,該方法對變量正態(tài)性不做要求,相比于Pearson相關(guān)其適用范圍更廣[12]。目前,Spearman秩相關(guān)分析在微生物組學(xué)領(lǐng)域中已被廣泛用于識別微生物與其他水平分子或臨床指標(biāo)之間的關(guān)聯(lián)性[13-15]。Spearman秩相關(guān)是一種非參數(shù)方法,適用范圍比Pearson更廣,但其中仍存在一些問題:(1)微生物數(shù)據(jù)中零觀測值出現(xiàn)的概率通常遠遠超出一般計數(shù)模型(如泊松回歸和負(fù)二項回歸模型)的預(yù)期范圍,如果繼續(xù)使用一般關(guān)聯(lián)分析方法可能得到錯誤結(jié)論。(2)Spearman秩相關(guān)方法作為一種非參數(shù)的方法,其統(tǒng)計檢驗效能較低。(3)只能識別具有單調(diào)性的相關(guān)關(guān)系,對于微生物關(guān)聯(lián)研究中的多種函數(shù)關(guān)系的非線性相關(guān)不敏感。

2.零膨脹模型

Lambert首次建立了零膨脹泊松回歸模型,它的基本思想是將計數(shù)數(shù)據(jù)中的零觀測值歸結(jié)于兩部分,一是來源于數(shù)據(jù)結(jié)構(gòu)的結(jié)構(gòu)零,二是由分布產(chǎn)生的抽樣零[16]。

零膨脹泊松回歸模型可以看作Bernoulli分布和Poisson分布組成的混合分布,其概率密度函數(shù)為:

Xu等通過大量的模擬實驗,比較了零膨脹模型與常用方法對具有零膨脹特征的數(shù)據(jù)進行建模的性能[17]。模擬研究表明,零膨脹模型可以較好控制一類錯誤率,并且具有更高的統(tǒng)計效能,對參數(shù)估計更準(zhǔn)確。

然而,泊松回歸模型假設(shè)事件發(fā)生的期望和方差相等,并且要求事件發(fā)生前后相互獨立,但微生物數(shù)據(jù)方差通常大于均數(shù),表現(xiàn)出過離散的現(xiàn)象,往往不符合這一假設(shè),這將導(dǎo)致模型參數(shù)估計值的標(biāo)準(zhǔn)誤偏小,參數(shù)Wald檢驗的假陽性率增加。為了解決零膨脹數(shù)據(jù)中的過度離散問題,Martin. Ridout研究了零膨脹負(fù)二項回歸模型,并證明了零膨脹負(fù)二項回歸比零膨脹泊松回歸模型更適合處理離散度高的數(shù)據(jù)[18]。零膨脹負(fù)二項回歸模型以負(fù)二項分布為基礎(chǔ),相比于零膨脹泊松回歸模型,能更好擬合不同離散度的數(shù)據(jù)。Wu等使用零膨脹負(fù)二項回歸模型發(fā)現(xiàn)表皮生長因子受體陽性的女性腸道菌群α多樣性較低,隨著乳腺癌惡化FirmicutesProteobacteria等微生物豐度增加[19]。Schwimmer等用同樣的方法研究了非酒精性脂肪肝患者的腸道微生物,發(fā)現(xiàn)炎癥相關(guān)菌群與非酒精性脂肪發(fā)病及其嚴(yán)重程度存在關(guān)聯(lián),Prevotella的豐度與肝纖維化有關(guān),并且使用菌群構(gòu)建了區(qū)分能力較好的判別模型[20]。

線性相關(guān)是容易解釋的相關(guān)形式,針對微生物數(shù)據(jù)特點,零膨脹模型對于微生物數(shù)據(jù)線性關(guān)系的檢驗效能顯著高于其他方法。但是,零膨脹模型在識別微生物和其他層面組學(xué)變量或一些臨床指標(biāo)間的復(fù)雜非線性關(guān)聯(lián)中存在局限性。

3.基于互信息的關(guān)聯(lián)分析方法

互信息 (mutual information,MI) 是來源于信息論中熵的概念,可用于度量兩個隨機變量之間依賴程度,表明一個隨機變量包含關(guān)于另一個隨機變量的信息量[21]?；バ畔⒌挠嬎阒?兩個隨機變量(X,Y)的聯(lián)合分布為p(x,y),邊緣分布為p(x),p(y),互信息I(X;Y)是聯(lián)合分布p(x,y)于邊緣分布p(x)p(y)的相對熵,即

根據(jù)熵的連鎖規(guī)則,有

H(X,Y)=H(X)+H(Y|X)=H(Y)+H(X|Y)

互信息對樣本的分布類型無特別要求,可以有效度量變量間線性關(guān)系和非線性關(guān)系,兩個隨機變量之間的互信息越大,則兩者之間的相關(guān)性就越強。傳統(tǒng)的互信息中要求變量是離散的且已知雙變量的聯(lián)合概率密度。微生物物種豐度是計數(shù)數(shù)據(jù),可以先對微生物進行分箱,再計算互信息。然而傳統(tǒng)互信息方法對于分箱的方式選擇具有敏感性,不同的分箱方式可能導(dǎo)致結(jié)果差異較大。

Kraskov提出使用基于K近鄰 (K-nearest neighbors,KNN) 的方法計算兩連續(xù)變量的互信息[22]。該方法無需知道概率密度函數(shù)形式,避免了對概率密度函數(shù)的估計,適用于非線性不規(guī)則分布的數(shù)據(jù)。其基本思想是在由隨機變量X和Y構(gòu)成的空間中首先找到給定樣本的k個近鄰樣本,再計算X和Y軸方向上距離小于K近鄰距離的樣本數(shù)目,并據(jù)此進行互信息估計。非參數(shù)方法除了K近鄰方法還包括核密度估計的互信息[23]。K近鄰和核密度估計方法將傳統(tǒng)的互信息適用范圍從離散變量擴展到連續(xù)變量。Jahagirdar等使用互信息方法構(gòu)建了代謝物-代謝物的關(guān)聯(lián)網(wǎng)絡(luò),并發(fā)現(xiàn)互信息方法與Pearson和Spearman方法分析的結(jié)果存在差異[24]。Numata等使用基于KNN的互信息方法分析擬南芥代謝物濃度數(shù)據(jù),發(fā)現(xiàn)互信息能夠檢測出Pearson相關(guān)系數(shù)無法發(fā)現(xiàn)的額外非線性相關(guān)[25]。K近鄰和核密度估計的互信息的缺點在于沒有合適的歸一化數(shù)據(jù)預(yù)處理方法,關(guān)聯(lián)強度大小不便于比較。

David N. Reshef在互信息方法的基礎(chǔ)上提出了最大信息系數(shù) (maximal information coefficient,MIC) 的方法[26]。MIC方法首先對變量X和Y構(gòu)成的散點圖進行網(wǎng)格化,并且求出不同網(wǎng)格化方案中的最大互信息值,再對最大的互信息值進行歸一化。MIC的計算方式如下:

Logares等將MIC方法用于分析海洋表層微生物,發(fā)現(xiàn)原核生物和海洋的溫度與氧氣有關(guān)[27]。Cao等提出了一種共表達網(wǎng)絡(luò)分析方法,即在加權(quán)基因共表達網(wǎng)絡(luò)分析 (weighted correlation network analysis,WGCNA) 的基礎(chǔ)上,結(jié)合了Pearson相關(guān)系數(shù)和最大信息系數(shù)作為配對基因之間的相似性度量,以此構(gòu)建共表達網(wǎng)絡(luò)。該方法能發(fā)現(xiàn)WGCNA方法所忽略的非線性相關(guān),識別更多潛在的包含癌癥信息的基因,具有更高的預(yù)測精度[28]。MIC的優(yōu)點是可用于各種函數(shù)形式的相關(guān)關(guān)系分析,在樣本量足夠大時能為不同類型單噪聲程度相似的相關(guān)關(guān)系給出相近系數(shù)。但MIC的缺點在于其統(tǒng)計效能較低,在進行大規(guī)模的探索分析時,需要較大的樣本量。

4.知識驅(qū)動的網(wǎng)絡(luò)構(gòu)建方法

Rob Knight提出了microbe-metabolite vectors (mmvec) 神經(jīng)元網(wǎng)絡(luò)方法,該方法專用于微生物組與代謝組數(shù)據(jù),通過學(xué)習(xí)代謝物和微生物共現(xiàn)概率,確定微生物與代謝物之間的關(guān)系[29]。Mmvec方法通過給定單個輸入微生物序列的情況下估計代謝物響應(yīng)強度,通過計算估計關(guān)聯(lián)與真實關(guān)聯(lián)的誤差對模型的權(quán)重進行調(diào)整,并進行內(nèi)部交叉驗證和迭代訓(xùn)練,預(yù)測微生物和代謝物的關(guān)聯(lián)關(guān)系。Mmvec相比于Pearson,Spearman和SparCC具有更高的F1得分,精確率和召回率。

Borenstein Lab提出了MIMOSA2模型,基于微生物在全基因組代謝模型(genome-scale metabolic models,GEMs)或KEGG中已知的代謝相關(guān)信息,計算微生物群落代謝能力 (community-wide metabolite potential,CMP)構(gòu)建代謝模型來估計群落組成對代謝物濃度的影響,并評估CMP與觀測到的代謝組特征的差異,從而發(fā)現(xiàn)相關(guān)的微生物和代謝物[30-32]。然而mmvec和MIMOSA2方法的適用范圍較小,僅能用于微生物和代謝物的關(guān)聯(lián)分析,從多組學(xué)的角度來看,微生物可能同其他多個水平的分子之間存在關(guān)聯(lián)甚至因果調(diào)控關(guān)系。

微生物組學(xué)關(guān)聯(lián)分析策略

多組學(xué)研究設(shè)計在成為探索疾病發(fā)生發(fā)展深層機制的有力手段的同時,也為統(tǒng)計分析帶來了新的機遇與挑戰(zhàn)。研究者們在進行微生物組學(xué)數(shù)據(jù)關(guān)聯(lián)分析過程中,常常忽略微生物組學(xué)數(shù)據(jù)零膨脹,過離散等數(shù)據(jù)特點,以及相應(yīng)統(tǒng)計分析方法原理與前提假設(shè),導(dǎo)致采用單一分析方法可能存在檢驗效能過低的問題,甚至出現(xiàn)錯誤結(jié)果。本文綜合上述分析方法特點及微生物組學(xué)數(shù)據(jù)特點,探索了微生物組學(xué)關(guān)聯(lián)分析策略(流程圖如圖1所示),具體總結(jié)如下:

圖1 微生物組學(xué)數(shù)據(jù)關(guān)聯(lián)分析策略流程圖

1.數(shù)據(jù)預(yù)處理。參數(shù)方法對數(shù)據(jù)分布有嚴(yán)格要求,根據(jù)使用的方法對數(shù)據(jù)進行相應(yīng)預(yù)處理是十分必要的。對于微生物成分?jǐn)?shù)據(jù),需要進行中心對數(shù)轉(zhuǎn)換:

其他水平分子和臨床指標(biāo)等一般的連續(xù)變量也建議進行處理(如標(biāo)準(zhǔn)化,對數(shù)比轉(zhuǎn)換,BOX-COX變換等),使變量分布為標(biāo)準(zhǔn)正態(tài)或近似標(biāo)準(zhǔn)正態(tài)分布。

2.數(shù)據(jù)特點描述:對原始的和預(yù)處理后的數(shù)據(jù)進行描述性統(tǒng)計分析,關(guān)注變量分布情況。此外,還應(yīng)描述微生物組中每個變量零觀測值比例,對零觀測值比例較高的變量,應(yīng)采用零膨脹模型和非參數(shù)的關(guān)聯(lián)分析方法。

3.線性關(guān)系識別。針對微生物數(shù)據(jù)特點,零膨脹模型對于微生物數(shù)據(jù)線性關(guān)系的檢驗效能顯著高于其他方法,并且一類錯誤率較低。同時,線性相關(guān)形式相較于非線性相關(guān)形式具有更好的生物可解釋性。使用零膨脹負(fù)二項回歸或零膨脹beta回歸模型(分別適用于微生物絕對豐度數(shù)據(jù)和相對豐度數(shù)據(jù))對數(shù)據(jù)進行模型擬合,根據(jù)回歸系數(shù)的大小以及回歸系數(shù)的顯著性得到兩個變量間的(偏)相關(guān)系數(shù)大小及統(tǒng)計檢驗顯著性,若相關(guān)系數(shù)檢驗的P值小于設(shè)定置信水平α,則認(rèn)為兩變量間存在線性相關(guān),否則按照存在非線性關(guān)聯(lián)進行進一步檢驗。

4.關(guān)聯(lián)分析。對上一步識別為線性相關(guān)的變量對,采用零膨脹模型估計相關(guān)關(guān)系強弱。對于不存在線性關(guān)系的變量對,可以通過非線性相關(guān)分析方法做進一步探索。MIC可以發(fā)現(xiàn)多種多樣的相關(guān)函數(shù)形式,Spearman相關(guān)分析方法雖然統(tǒng)計檢驗效能高于MIC,但只能發(fā)現(xiàn)存在單調(diào)關(guān)系的變量。因此可同時采用MIC和Spearman的方法分別進行非線性相關(guān)關(guān)系識別,計算變量間的最大信息系數(shù)和相關(guān)系數(shù),并給出系數(shù)檢驗的P值,若P值小于置信水平α,則認(rèn)為兩變量間存在非線性相關(guān)。最后,給出線性和非線性相關(guān)的R和P值。

5.多重假設(shè)檢驗問題處理。高維情境常涉及多重檢驗問題,傳統(tǒng)的統(tǒng)計檢驗方法中會產(chǎn)生大量假陽性結(jié)果,大大增加假陽性錯誤[33]。對于多重檢驗,在檢驗水平固定的情況下,隨著檢驗次數(shù)的增加,至少犯一次假陽性錯誤的概率將會趨于1。因此需要考慮對假設(shè)檢驗結(jié)果的P值進行校正,將校正后的P值控制在某一固定水平及以下。目前常用校正方法包括Bonferroni校正法及FDR校正法[34]。對于第二步和第三步過程中的多重假設(shè)檢驗問題,采用FDR的方法校正P值的閾值。

由于微生物組學(xué)數(shù)據(jù)的復(fù)雜特性,上述分析手段能在一定程度上解決組學(xué)數(shù)據(jù)統(tǒng)計分析問題,但仍存在局限性。由于各變量間存在多種非線性關(guān)系,實際分析中可能需要結(jié)合多種網(wǎng)絡(luò)構(gòu)建方法[35]?；诨バ畔⒌姆椒m然適用于多種非線性關(guān)系并對數(shù)據(jù)分布沒有要求,但是在小樣本條件下檢驗效能較低,如何提高關(guān)聯(lián)分析方法適用范圍和檢驗效能需要進一步研究。此外,橫斷面研究對于探索微生物與其他水平分子的調(diào)控及因果關(guān)系可能存在證據(jù)不充分的問題,縱向研究中的研究策略和研究方法需要進一步探討。以上三個關(guān)鍵科學(xué)問題的解決將會對微生物組學(xué)關(guān)聯(lián)分析提供新的思路與契機。