微生物組學關聯分析方法和策略*

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

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

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

微生物組學數據特點

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

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

表1 微生物組學數據分布

微生物組學關聯分析常用統計分析方法

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

表2 關聯分析方法原理及優點

1.傳統關聯分析方法及其擴展方法

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

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

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

2.零膨脹模型

Lambert首次建立了零膨脹泊松回歸模型,它的基本思想是將計數數據中的零觀測值歸結于兩部分,一是來源于數據結構的結構零,二是由分布產生的抽樣零[16]。

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

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

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

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

3.基于互信息的關聯分析方法

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

根據熵的連鎖規則,有

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

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

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

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

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

4.知識驅動的網絡構建方法

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

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

微生物組學關聯分析策略

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

圖1 微生物組學數據關聯分析策略流程圖

1.數據預處理。參數方法對數據分布有嚴格要求,根據使用的方法對數據進行相應預處理是十分必要的。對于微生物成分數據,需要進行中心對數轉換:

其他水平分子和臨床指標等一般的連續變量也建議進行處理(如標準化,對數比轉換,BOX-COX變換等),使變量分布為標準正態或近似標準正態分布。

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

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

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

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

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