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

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

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

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

原理與方法

1.條件高斯貝葉斯網(wǎng)絡(luò)

貝葉斯網(wǎng)絡(luò)是由代表隨機(jī)變量的節(jié)點(diǎn)和變量間關(guān)系的邊構(gòu)成的有向無環(huán)圖。對(duì)于兼有離散變量和連續(xù)變量的數(shù)據(jù)，如本研究實(shí)例數(shù)據(jù)包含組學(xué)數(shù)據(jù)和分類結(jié)局，可以構(gòu)建條件高斯貝葉斯網(wǎng)絡(luò)。條件高斯貝葉斯網(wǎng)絡(luò)規(guī)定離散變量的父節(jié)點(diǎn)不能為連續(xù)變量，通過分別對(duì)連續(xù)型變量和離散型變量計(jì)算概率分布，構(gòu)成所有隨機(jī)變量的聯(lián)合概率分布，如下式所示：

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

(1)

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

離散節(jié)點(diǎn)的似然函數(shù)為：

(2)

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

連續(xù)節(jié)點(diǎn)的似然函數(shù)為：

(3)

其中zi是變量yi的連續(xù)型父節(jié)點(diǎn)的值，βij是當(dāng)變量yi的離散型父節(jié)點(diǎn)分類為j時(shí)zi變量對(duì)yi變量回歸的參數(shù)向量，而τij為精度，是方差的倒數(shù)。給定數(shù)據(jù)集后，結(jié)合似然函數(shù)與先驗(yàn)分布估計(jì)不同網(wǎng)絡(luò)結(jié)構(gòu)的后驗(yàn)概率，即所有變量聯(lián)合分布的后驗(yàn)概率，后驗(yàn)概率取對(duì)數(shù)后即為貝葉斯狄利克雷評(píng)分(Bayesian Dirichlet score，BD)。

首先基于貝葉斯因子進(jìn)行變量粗篩。貝葉斯因子是某個(gè)自變量和因變量相關(guān)聯(lián)與二者獨(dú)立時(shí)的對(duì)數(shù)似然值之差，貝葉斯因子閾值大于1即表示有中等證據(jù)支持節(jié)點(diǎn)之間的關(guān)聯(lián)，大于3即代表有較強(qiáng)的證據(jù)支持關(guān)聯(lián)[10]。通過設(shè)置貝葉斯因子閾值初步篩選可能對(duì)分類造成影響的變量，達(dá)到降維目的。

貝葉斯網(wǎng)絡(luò)的結(jié)構(gòu)學(xué)習(xí)采用基于評(píng)分的搜索算法，即通過BD評(píng)分度量網(wǎng)絡(luò)結(jié)構(gòu)與樣本數(shù)據(jù)的擬合程度;然后采用K2算法來尋找評(píng)分最高即最優(yōu)的網(wǎng)絡(luò)結(jié)構(gòu)。K2算法是根據(jù)變量的似然值對(duì)節(jié)點(diǎn)進(jìn)行排序，對(duì)每個(gè)變量根據(jù)序號(hào)從前到后選擇父節(jié)點(diǎn)，從而減少節(jié)點(diǎn)選擇的可能組合，降低計(jì)算負(fù)擔(dān)。貝葉斯網(wǎng)絡(luò)的參數(shù)學(xué)習(xí)是基于獲得的網(wǎng)絡(luò)結(jié)構(gòu)，從訓(xùn)練樣本獲得各節(jié)點(diǎn)的條件概率分布表，采用極大似然法估計(jì)網(wǎng)絡(luò)參數(shù)。分類判別是計(jì)算分類變量的邊際概率估計(jì)完成的。

2.模擬研究

模擬數(shù)據(jù)分別考慮不同相關(guān)程度和不同稀疏水平，同時(shí)考慮線性和非線性關(guān)系。每類數(shù)據(jù)均模擬200個(gè)連續(xù)型自變量以及一個(gè)二分類因變量。

(1)不同相關(guān)系數(shù)的數(shù)據(jù)生成：首先生成標(biāo)準(zhǔn)正態(tài)分布變量X1，然后通過回歸公式Xn=βX1+ε生成與X1相關(guān)的99個(gè)變量，再生成100個(gè)標(biāo)準(zhǔn)正態(tài)分布變量，合并后即為具有一定相關(guān)性的模板數(shù)據(jù)，以該數(shù)據(jù)集的協(xié)方差矩陣，以及0為均數(shù)，生成模擬數(shù)據(jù)。不同模擬數(shù)據(jù)的變量間相關(guān)性高低通過回歸公式的β來控制，本研究?jī)深悢?shù)據(jù)集的β值為0.5和2。

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

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

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

3.實(shí)例分析

(1)數(shù)據(jù)描述

本研究所使用數(shù)據(jù)來自于美國國家代謝組學(xué)數(shù)據(jù)存儲(chǔ)庫NMDR (national metabolomics data repository)網(wǎng)站，網(wǎng)址https://www.metabolomicsworkbench.org，項(xiàng)目號(hào)為PR000300。該數(shù)據(jù)可以通過項(xiàng)目DOI：10.21228/M88C86直接獲得。該數(shù)據(jù)庫由NIH，U2C- DK119886資助?？梢酝ㄟ^以下網(wǎng)址進(jìn)入:https://www.metabolomicsworkbench.org/about/howtocite.php.

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

(2)實(shí)例數(shù)據(jù)分析方法

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

采用有放回重抽樣擴(kuò)增樣本例數(shù)少的一組，獲得平衡數(shù)據(jù)。根據(jù)最佳模型參數(shù)，對(duì)平衡數(shù)據(jù)計(jì)算兩個(gè)方法的最終AUC，并進(jìn)行permutation檢驗(yàn)[12]。

結(jié) 果

1.模擬研究結(jié)果

(1)不同相關(guān)系數(shù)模擬數(shù)據(jù)

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

圖1 不同相關(guān)程度模擬數(shù)據(jù)分類評(píng)價(jià)結(jié)果

表1 不同相關(guān)程度模擬數(shù)據(jù)分類AUC

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

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

表2 不同相關(guān)程度模擬數(shù)據(jù)分類AUC

圖2 不同稀疏水平模擬數(shù)據(jù)分類評(píng)價(jià)結(jié)果

2.實(shí)例分析結(jié)果

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

討 論

本研究基于模擬數(shù)據(jù)和實(shí)例數(shù)據(jù)比較了PLSDA和CGBN的分類判別能力。研究發(fā)現(xiàn)，對(duì)于獲取成本較高的小樣本高維度代謝組學(xué)數(shù)據(jù)，并且在代謝物之間非線性相關(guān)時(shí)CGBN可能是值得探討的新方法。

模擬結(jié)果顯示在中小樣本量時(shí)，無論變量之間相關(guān)系數(shù)或者稀疏程度的高低，CGBN均是優(yōu)于PLSDA的方法。而樣本量較大時(shí)，對(duì)于自變量與因變量線性相關(guān)的數(shù)據(jù)或者低相關(guān)非線性數(shù)據(jù)，CGBN仍然能取得更好的分類效果。

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

本研究采用預(yù)先確定的貝葉斯因子閾值來對(duì)節(jié)點(diǎn)進(jìn)行篩選，理論上可以調(diào)節(jié)進(jìn)入模型的節(jié)點(diǎn)數(shù)，建立較簡(jiǎn)約的模型，一定程度避免過擬合的發(fā)生[13]。本文模擬研究時(shí)不同樣本量情形下均設(shè)貝葉斯因子閾值為1，這是一個(gè)較寬松的標(biāo)準(zhǔn)，因此樣本量較大時(shí)，可能會(huì)不能有效篩選變量，模型不夠簡(jiǎn)約，導(dǎo)致過擬合而影響模型外推性。因此，實(shí)際應(yīng)用中需根據(jù)具體的研究目的和數(shù)據(jù)維度來設(shè)定合適的貝葉斯因子閾值[14]。

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

模擬研究發(fā)現(xiàn)CGBN方法的計(jì)算負(fù)擔(dān)相對(duì)較大，因此，對(duì)于研究者可以結(jié)合具體的研究目標(biāo)和擁有的計(jì)算資源決定方法的選擇，在不考慮計(jì)算負(fù)擔(dān)而追求更高的分類準(zhǔn)確率時(shí)建議使用CGBN。

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