基于條件高斯貝葉斯網絡的代謝組學數據分類預測研究*

戴品遠 余小金△ 謝緯華 趙 超 劉 冉 尹立紅 陳炳為

【提 要】 目的 探討條件高斯貝葉斯網絡(conditional Gaussian Bayesian network，CGBN)在代謝組學數據的分類判別中的應用。方法 通過模擬研究與實際代謝組學數據分析，比較CGBN與偏最小二乘判別分析(partial least squares discriminant analysis，PLSDA)在不同相關程度和不同稀疏水平的高維數據及線性相關與非線性等情形時的分類判別性能，評價指標采用ROC曲線下面積(area under curve，AUC)和平均計算時間。結果 模擬研究結果表明，變量之間低相關且樣本量不大于200時CGBN分類判別AUC高于PLSDA。在自變量與因變量非線性相關且小樣本情況下CGBN分類判別AUC同樣高于PLSDA。實例數據分析結果顯示CGBN和PLSDA分類判別的AUC分別為0.997，0.975。CGBN的計算時間要遠高于PLSDA。結論 在不受計算負擔限制的情形下，CGBN是代謝組學數據典型分析方法的一種可行的替代方法，值得進一步研究。

近年來代謝組學(metabolomics)數據的分析方法迅速應用于包括疾病診斷，療效評價等領域[1-7]。目前常用的代謝組學數據分類方法主要包括偏最小二乘判別分析(partial least squares discriminant analysis，PLSDA)、支持向量機(support vector machine，SVM)、神經網絡(artificial neural networks，ANNs)等。而基于貝葉斯網絡的分析方法在代謝組學數據的分析中越來越受研究者關注。與應用最為普遍的PLSDA相比，貝葉斯網絡分析可識別代謝物和表型之間的非線性關系，并結合貝葉斯先驗概率，減輕過擬合程度；識別有差異的代謝物子網[8]，因此，成為目前應用領域的研究熱點。

Kelly[9]等基于實際數據分析發現，即使在無信息先驗情況下，條件高斯貝葉斯網絡(conditional Gaussian Bayesian network，CGBN)仍優于經典偏最小二乘方法判別分析，并且具有更小的過擬合風險。但對于更一般的數據特征，基于貝葉斯網絡思維的分類判別的性能仍需要深入探討。本研究基于模擬數據探討CGBN與PLSDA相比較的分類預測能力和計算效能，并將兩種方法應用于肥胖婦女人群的2型糖尿病病例與對照分類判別中。

原理與方法

1.條件高斯貝葉斯網絡

貝葉斯網絡是由代表隨機變量的節點和變量間關系的邊構成的有向無環圖。對于兼有離散變量和連續變量的數據，如本研究實例數據包含組學數據和分類結局，可以構建條件高斯貝葉斯網絡。條件高斯貝葉斯網絡規定離散變量的父節點不能為連續變量，通過分別對連續型變量和離散型變量計算概率分布，構成所有隨機變量的聯合概率分布，如下式所示：

p(Δ)f(Ψ|Δ)=∏x∈Δp(x|π(x))∏y∈Ψf(y|π(y))

(1)

其中Δ為離散變量，Ψ為連續變量，p(x|π(x))為給定變量x的父節點π(x)后的條件概率，而f(y|π(y))表示給定變量y的父節點π(y)后的條件分布。

離散節點的似然函數為：

(2)

其中nijk表示離散變量xi分類為k，父節點π(xi)為分類為j時的樣本數量，而αijk是離散變量xi分類為k，父節點π(xi)為分類為j時狄利克雷分布中先驗頻數的超參數。

連續節點的似然函數為：

(3)

其中zi是變量yi的連續型父節點的值，βij是當變量yi的離散型父節點分類為j時zi變量對yi變量回歸的參數向量，而τij為精度，是方差的倒數。給定數據集后，結合似然函數與先驗分布估計不同網絡結構的后驗概率，即所有變量聯合分布的后驗概率，后驗概率取對數后即為貝葉斯狄利克雷評分(Bayesian Dirichlet score，BD)。

首先基于貝葉斯因子進行變量粗篩。貝葉斯因子是某個自變量和因變量相關聯與二者獨立時的對數似然值之差，貝葉斯因子閾值大于1即表示有中等證據支持節點之間的關聯，大于3即代表有較強的證據支持關聯[10]。通過設置貝葉斯因子閾值初步篩選可能對分類造成影響的變量，達到降維目的。

貝葉斯網絡的結構學習采用基于評分的搜索算法，即通過BD評分度量網絡結構與樣本數據的擬合程度;然后采用K2算法來尋找評分最高即最優的網絡結構。K2算法是根據變量的似然值對節點進行排序，對每個變量根據序號從前到后選擇父節點，從而減少節點選擇的可能組合，降低計算負擔。貝葉斯網絡的參數學習是基于獲得的網絡結構，從訓練樣本獲得各節點的條件概率分布表，采用極大似然法估計網絡參數。分類判別是計算分類變量的邊際概率估計完成的。

2.模擬研究

模擬數據分別考慮不同相關程度和不同稀疏水平，同時考慮線性和非線性關系。每類數據均模擬200個連續型自變量以及一個二分類因變量。

(1)不同相關系數的數據生成：首先生成標準正態分布變量X1，然后通過回歸公式Xn=βX1+ε生成與X1相關的99個變量，再生成100個標準正態分布變量，合并后即為具有一定相關性的模板數據，以該數據集的協方差矩陣，以及0為均數，生成模擬數據。不同模擬數據的變量間相關性高低通過回歸公式的β來控制，本研究兩類數據集的β值為0.5和2。

(2)不同稀疏水平的數據生成：以稀疏性P控制變量之間的稀疏水平，從而獲得數據的結構矩陣，然后通過這個結構矩陣按照自由度為9，Σ為單位矩陣的G-Wishart分布生成相應的協方差矩陣，最后以0為均數生成模擬數據。本研究所使用的P為0.3和0.8。

每種類型都生成樣本量為50，100，200和500的模擬數據，再對數據的變量進行線性或者非線性的變換，得到變換后的變量hi(Xi)及其組合g(Xi)，進而得到每個樣本的分類因變量。以1：1的比例隨機分為訓練集和測試集。

使用PLSDA和CGBN對模擬數據進行分類判別；PLSDA分析選擇因變量方差解釋累計貢獻率達90%的成分進行分析，而貝葉斯網絡分析時將貝葉斯因子閾值設為1對變量進行初篩。每次判別分析后記錄兩個模型的運算時間，并建立混淆矩陣，計算受試者工作特征曲線下面積，AUC的值。對每個樣本量均重復1000次模擬實驗，計算AUC和運算時間的均數、標準差，獲得AUC的置信區間。采用MATLAB(R2019a)CGBayesNets包學習CGBN并進行分類判別，結構學習使用的離散變量先驗超參數α=10，連續變量先驗方差∑2=1，限制最大父節點數為 2。本研究使用R的mixOmics包來實現PLSDA。全部運算均在2.40 GHz Intel Xeon CPU計算機上進行。

3.實例分析

(1)數據描述

本研究所使用數據來自于美國國家代謝組學數據存儲庫NMDR (national metabolomics data repository)網站，網址https://www.metabolomicsworkbench.org，項目號為PR000300。該數據可以通過項目DOI：10.21228/M88C86直接獲得。該數據庫由NIH，U2C- DK119886資助。可以通過以下網址進入:https://www.metabolomicsworkbench.org/about/howtocite.php.

本研究所用代謝組學公共數據庫研究對象為患2型糖尿病(n=44)和非糖尿病(n=12)的非裔美國體重超標婦女。該樣本年齡范圍在19.3到87.1歲之間，基于GC-TOF質譜法分析了血漿樣本的代謝組學信息后，使用BinBase數據庫過濾，共360種代謝物通過了嚴格的質量控制措施。

(2)實例數據分析方法

首先使用PLSDA通過5折交叉驗證，分別計算1到20個成分時模型的AUC值，以因變量方差解釋累計貢獻率達到90%及以上和獲得最大AUC確定模型參數[11]；貝葉斯網絡分析首先通過1～20的貝葉斯因子閾值篩選不同的變量，每個貝葉斯因子閾值篩選后的數據采用重復數為100的bootstrap法分別學習CGBN，基于預先給定的頻率閾值(本文采用0.3～0.8等6個閾值)，以AUC最大原則選擇最佳模型參數，確定最后的貝葉斯網絡結構模型。

采用有放回重抽樣擴增樣本例數少的一組，獲得平衡數據。根據最佳模型參數，對平衡數據計算兩個方法的最終AUC，并進行permutation檢驗[12]。

結 果

1.模擬研究結果

(1)不同相關系數模擬數據

對模擬數據的設定特征進行檢查，低相關數據集存在相關的自變量之間的平均相關系數為0.171(0.170，0.172)，高相關數據集的平均相關系數為0.801 (0.799，0.802)。模擬研究結果如表1和圖1所示。在變量間低相關的情況下，樣本量小于等于200時CGBN的AUC要高于PLSDA。在高相關非線性時，兩種方法較為接近，小樣本情況下CGBN的AUC較高，而高相關線性時CGBN的AUC較高。經計算，PLSDA的平均計算時間為0.0109(0.0108，0.0111)秒，而CGBN的平均計算時間為4.1800(4.1661，4.1937)秒。

圖1 不同相關程度模擬數據分類評價結果

表1 不同相關程度模擬數據分類AUC

(2)不同稀疏水平模擬數據

低稀疏數據集的稀疏性為0.804，高稀疏數據集的稀疏性為0.303，稀疏性為數據的精度矩陣中不為0的元素所占比例，數值越小則越稀疏。模擬研究結果如表2和圖2所示，當自變量與因變量是非線性相關時，若樣本量小于等于100，則CGBN的AUC要高于PLSDA，而對于線性相關的數據，則在較小樣本時兩種方法AUC接近，隨樣本量增加PLSDA的判別能力優于CGBN。PLSDA的平均計算時間為0.0294(0.0292，0.0296)秒，而CGBN的平均計算時間為1.9287(1.9149，1.9424)秒。

表2 不同相關程度模擬數據分類AUC

圖2 不同稀疏水平模擬數據分類評價結果

2.實例分析結果

以4為貝葉斯因子閾值篩選變量，使用CGBN對平衡數據分類的AUC為0.997，靈敏度為0.977，特異度為1，誤判率為0.011。選擇兩個成分進行PLSDA，AUC為0.975，靈敏度為0.955，特異度為0.977，誤判率為0.034。permutation檢驗表明，PLSDA與CGBN的AUC計算得到的P值均小于0.001。

討 論

本研究基于模擬數據和實例數據比較了PLSDA和CGBN的分類判別能力。研究發現，對于獲取成本較高的小樣本高維度代謝組學數據，并且在代謝物之間非線性相關時CGBN可能是值得探討的新方法。

模擬結果顯示在中小樣本量時，無論變量之間相關系數或者稀疏程度的高低，CGBN均是優于PLSDA的方法。而樣本量較大時，對于自變量與因變量線性相關的數據或者低相關非線性數據，CGBN仍然能取得更好的分類效果。

本研究實例數據為通過體質指數和年齡進行匹配所獲得的單性別病例對照樣本，生活在相對較小的地理空間內，遺傳變異較低，并且具有共同的飲食攝入模式，人群生物代謝信號噪聲較低，適合進行代謝組學研究。分析結果顯示CGBN的AUC稍高于PLSDA，且經permutation檢驗不存在過擬合風險。因此，本研究對實例數據的分析也顯示CGBN方法具有分類準確性優勢。但由于樣本量較小且不平衡，結果解釋需要更豐富的數據支持。

本研究采用預先確定的貝葉斯因子閾值來對節點進行篩選，理論上可以調節進入模型的節點數，建立較簡約的模型，一定程度避免過擬合的發生[13]。本文模擬研究時不同樣本量情形下均設貝葉斯因子閾值為1，這是一個較寬松的標準，因此樣本量較大時，可能會不能有效篩選變量，模型不夠簡約，導致過擬合而影響模型外推性。因此，實際應用中需根據具體的研究目的和數據維度來設定合適的貝葉斯因子閾值[14]。

本研究分別模擬不同的相關程度和不同的稀疏水平數據，沒有探討兩參數交叉條件下的方法性能。在生成不同稀疏程度的模擬數據時得到的相關系數取值較為分散，因此，本研究所展示的稀疏水平對分類效果的影響可能受不同數據集的相關系數的影響。另一方面，本研究模擬分析時將樣本隨機分為訓練集和驗證集，在小樣本情形下，與交叉驗證相比較可能無法準確描述模型的外推能力[15]，但可以避免大樣本時模擬運行的計算負擔。因此，可以考慮在不同的模型驗證方案下進一步探討模型的實際外推性能。

模擬研究發現CGBN方法的計算負擔相對較大，因此，對于研究者可以結合具體的研究目標和擁有的計算資源決定方法的選擇，在不考慮計算負擔而追求更高的分類準確率時建議使用CGBN。

本研究采用基于評分的搜索算法完成CGBN結構學習。但是，在高維低樣本量的情況下可能會出現局部最優的情況[16]，因此需要進行算法改進。一些學者提出可以將基于評分的算法與基于條件獨立的算法相結合，比如最大最小爬山算法(max-min hill-climbing，MMHC)[17]，CB算法[18]等。此外，基于圖論的分析思維可能為代謝組學的CGBN圖模型的全局分析和局部分析提供更多信息。例如，可以分析與疾病發生風險有關系的一簇代謝物間的作用模式，即代謝子網，而不僅僅是單個代謝物的定量增減或者尋找代謝網絡的關鍵變量等，都將是未來的研究方向。