第三軍醫(yī)大學(xué)軍事預(yù)防醫(yī)學(xué)院衛(wèi)生統(tǒng)計(jì)學(xué)教研室（400038）

吳小姣 李高明 易大莉 劉 嶺 張彥琦 易 東 伍亞舟△

·方法介紹·

基因表達(dá)譜的非參缺失森林填補(bǔ)算法研究*

第三軍醫(yī)大學(xué)軍事預(yù)防醫(yī)學(xué)院衛(wèi)生統(tǒng)計(jì)學(xué)教研室（400038）

吳小姣 李高明 易大莉 劉 嶺 張彥琦 易 東 伍亞舟△

目的評(píng)估不同基因表達(dá)譜數(shù)據(jù)集下，多種算法在缺失數(shù)據(jù)中的填補(bǔ)效果，并分析其對(duì)無監(jiān)督基因表達(dá)譜聚類的影響。方法在不同缺失比例的情況下，分別采用非參缺失森林填補(bǔ)法、貝葉斯線性回歸法、蒙特卡洛多重填補(bǔ)法和K鄰近填補(bǔ)法進(jìn)行填補(bǔ)，通過均方根誤差（NRMSE）和聚類準(zhǔn)確率（F值）兩個(gè)指標(biāo)來評(píng)估不同方法的填補(bǔ)效能和聚類效果，并用模擬數(shù)據(jù)集進(jìn)行測試和乳腺癌數(shù)據(jù)集進(jìn)行驗(yàn)證。結(jié)果隨著缺失比例的增加，四種填補(bǔ)方法的NRMSE都逐漸上升；任意缺失比例下，相比于其他三種方法非參缺失森林填補(bǔ)法的填補(bǔ)優(yōu)勢明顯。缺失比例為5%、10%、20%和30%的乳腺癌數(shù)據(jù)集，非參缺失森林填補(bǔ)法的NRMSE依次為0．1951（95%CI，0．1945～0．1953）、0．2776（95%CI，0．2783～2791）、0．4003（95%CI，0．3986～0．4002）和0．4974（95%CI，0．4658～0．5104）；聚類效果的準(zhǔn)確率為1．0、0．91、0．88和0．82。結(jié)論非參缺失森林填補(bǔ)算法實(shí)現(xiàn)簡單，對(duì)數(shù)據(jù)集的要求較低，比傳統(tǒng)填補(bǔ)算法具有更好的穩(wěn)定性和精確度，可保留較多的基因信息供后續(xù)的功能聚類等分析。

基因表達(dá)譜 缺失數(shù)據(jù) 缺失森林法 聚類

生物醫(yī)學(xué)研究中基因表達(dá)微陣列是一種強(qiáng)有力的工具，但現(xiàn)存的很多分析方法都要求微陣列的數(shù)據(jù)是完整的。由于存在多種原因，如不充分的實(shí)驗(yàn)方案，圖像損壞，芯片上的灰塵或劃痕等，使得實(shí)際上獲得的數(shù)據(jù)陣列通常是有缺失的，這在一定程度上影響了數(shù)據(jù)后續(xù)分析結(jié)果的準(zhǔn)確性和可靠性，如差異表達(dá)基因的篩選、基因功能聚類、基因調(diào)控網(wǎng)絡(luò)建立和生物標(biāo)志物檢測等。目前芯片缺失數(shù)據(jù)填補(bǔ)估計(jì)方法的文獻(xiàn)較多，可以大致分為四類：（1）局部算法：K鄰近距離法［1］、局部最小二乘法［2］等；（2）全局算法：奇異值分解法［3］、貝葉斯填補(bǔ)算法［4］等；（3）混合算法：linC-mb［5］；（4）利用生物信息輔助算法：POCS［6］、HAI填補(bǔ)［7］等。這些填補(bǔ)方法大多屬于參數(shù)統(tǒng)計(jì)方法，且都要求數(shù)據(jù)集的分布已知。實(shí)際上，基因表達(dá)譜數(shù)據(jù)集往往具有復(fù)雜的數(shù)據(jù)結(jié)構(gòu)且無任何先驗(yàn)知識(shí)，非參數(shù)模型方法對(duì)此卻能取得很好的效果；同時(shí)針對(duì)不同數(shù)據(jù)集的不同分析目的，將多種方法同時(shí)進(jìn)行比較的文獻(xiàn)較少，其研究尚有較大空間。本文介紹的非參缺失森林填補(bǔ)算法即為一種非參數(shù)統(tǒng)計(jì)方法，首次將其應(yīng)用于基因表達(dá)譜缺失數(shù)據(jù)的填補(bǔ)，并將其與常用的幾種填補(bǔ)方法（如貝葉斯線性回歸法［4］、蒙特卡洛多重填補(bǔ)法［8］和K鄰近填補(bǔ)法［2］）的填補(bǔ)效果進(jìn)行比較，最后分析各種填補(bǔ)方法對(duì)無監(jiān)督基因表達(dá)譜聚類的影響，為同類研究提供方法學(xué)借鑒。

理論與算法

1．非參缺失森林的填補(bǔ)方法

隨機(jī)森林算法［9］要求應(yīng)變量是完整的，才能訓(xùn)練出森林，Stekhoven在此基礎(chǔ)上進(jìn)行改進(jìn)，提出了缺失森林算法［10］，它可以直接用已觀測到的完整部分?jǐn)?shù)據(jù)集訓(xùn)練出的隨機(jī)森林來預(yù)測缺失值，而不依賴于應(yīng)變量的完整性。

假定數(shù)據(jù)集X＝（X1，X2，…，Xp）是一個(gè)N×P維的矩陣（N個(gè)基因，P個(gè)樣本），將其中任意一個(gè)可能含缺失數(shù)據(jù)的變量記為XS。應(yīng)變量和自變量的觀測值、缺失數(shù)據(jù)分別記為Yobs、Ymis和Xobs、Xmis。

具體的填補(bǔ)步驟如下：首先，用均數(shù)或其他填補(bǔ)方法對(duì)X的所有缺失值作初步的猜測，并將變量XS按缺失值的數(shù)量升序排列，令這個(gè)初步填補(bǔ)后的矩陣為Xold。對(duì)每一個(gè)變量XS，缺失森林算法的填補(bǔ)過程為：

（1）首先用應(yīng)變量Yobs和自變量Xobs擬合一個(gè)隨機(jī)森林；

（2）然后將Xmis作為特征變量輸入，用訓(xùn)練后的隨機(jī)森林來預(yù)測缺失數(shù)據(jù)Ymis，令新預(yù)測填補(bǔ)后得到的矩陣為Xnew；

（3）重復(fù)此填補(bǔ)過程，直到符合停止標(biāo)準(zhǔn)γ，即新填補(bǔ)的數(shù)據(jù)矩陣Xnew和前一個(gè)數(shù)據(jù)矩陣Xold的差值首次開始增加時(shí)；連續(xù)變量N間的差值定義為

2．缺失填補(bǔ)的效果評(píng)價(jià)

任何一種填補(bǔ)方法，都有各自的優(yōu)點(diǎn)和缺點(diǎn)。通常情況下，采用標(biāo)準(zhǔn)化均方根誤差（normalized root mean square error，NRMSE）作為不同方法填補(bǔ)效果的評(píng)價(jià)指標(biāo)［11］：

其中Ri為估計(jì)值，Ii為原始值，Std（Ii）為原始值的標(biāo)準(zhǔn)差。NRMSE的值越小表示其估計(jì)的越準(zhǔn)確、性能越好，反之結(jié)果越差。

3．聚類的效能評(píng)價(jià)

用層次聚類法對(duì)四種方法填補(bǔ)的完整數(shù)據(jù)集進(jìn)行聚類，并對(duì)聚類結(jié)果進(jìn)行評(píng)價(jià)。層次聚類法產(chǎn)生一個(gè)嵌套聚類的層次，算法最多包含N步，在第t步執(zhí)行的操作就是在前t－1步的聚類基礎(chǔ)上生成新聚類。假定對(duì)N個(gè)對(duì)象進(jìn)行聚類，層次聚類法［12－14］過程如下：

（1）初始時(shí)共有N類，每個(gè)類有一個(gè)對(duì)象構(gòu)成。令序號(hào)m＝0，l（m）＝0。

（2）在D中尋找最小距離d［r，s］＝min d［（i），（j）］。

（3）將兩個(gè)類r和類s合并成一個(gè)新類（R，S），令m＝m+1，L（m）＝d［r，s］。

（4）更新距離矩陣D：將表示類r和類s的行列刪除，同時(shí)加入表示新類（r，s）的行列；同時(shí)定義新類（r，s）與舊類（k）的距離為d［（k），（r，s）］＝min（d［（k），（r）］，d［（k），（s）］）。

（5）重復(fù)（2）～（4）步，直到所有對(duì)象合并成一個(gè)類為止。

在聚類的過程中，每次抽取一個(gè)填補(bǔ)缺失數(shù)據(jù)的方法為檢驗(yàn)樣本，以完整數(shù)據(jù)集的樣本構(gòu)成訓(xùn)練集，用訓(xùn)練集訓(xùn)練分類器，然后對(duì)檢驗(yàn)樣本進(jìn)行檢驗(yàn)，分別記錄下每個(gè)支持向量機(jī)在檢驗(yàn)樣本陽性類和陰性類的真陽性數(shù)（TP），真陰性數(shù)（TN），假陽性數(shù)（FP），假陰性數(shù)（FN）。一般用F值［15］方法對(duì)聚類的效果進(jìn)行評(píng)價(jià)：

其中，P＝TP/（TP+FP）；R＝TP/（TP+FN）；β為偽錯(cuò)誤的概率。F值越大表示其聚類效能越好，反之越差。

數(shù)據(jù)集

1．?dāng)?shù)據(jù)來源

本實(shí)驗(yàn)采用兩個(gè)基因表達(dá)數(shù)據(jù)集，第一個(gè)使用R軟件的ARTIVA包模擬一個(gè)多元正態(tài)分布的表達(dá)譜數(shù)據(jù)集，表示1024個(gè)基因在15個(gè)實(shí)驗(yàn)水平下的不同表達(dá)。第二個(gè)數(shù)據(jù)集來自GEO數(shù)據(jù)庫上公開發(fā)表的乳腺癌基因表達(dá)譜數(shù)據(jù)［16］，該數(shù)據(jù)集為6365個(gè)基因，15個(gè)實(shí)驗(yàn)樣本，兩個(gè)數(shù)據(jù)集都為非時(shí)間序列型結(jié)構(gòu)。

2．統(tǒng)計(jì)分析

分別對(duì)模擬和乳腺癌數(shù)據(jù)集，采用統(tǒng)計(jì)軟件包R3．2．4編程，按照一定百分比（如5%、10%、20%、30%）產(chǎn)生隨機(jī)性缺失數(shù)據(jù)，在統(tǒng)計(jì)軟件R下分別使用非參缺失森林法、貝葉斯線性回歸法、蒙特卡洛多重填補(bǔ)法和K鄰近法對(duì)缺失的乳腺癌表達(dá)譜數(shù)據(jù)集進(jìn)行填補(bǔ)，并進(jìn)行基因功能聚類分析的效果評(píng)估。需要加載的程序包有：affy、compositions、mice、missForest、impute、hclust、cutree。

結(jié) 果

1．基于均方根的填補(bǔ)效果評(píng)價(jià)

四種算法的填補(bǔ)效果如圖1所示。無論使用哪種填補(bǔ)方法，NRMSE的值都會(huì)隨著缺失比例的增加而逐漸上升。如乳腺癌數(shù)據(jù)集在缺失比例為10%時(shí)，非參缺失森林法、蒙特卡洛多重填補(bǔ)法、K鄰近填補(bǔ)法和貝葉斯線性回歸法的NRMSE依次為0．2671、0．3202、0．3190和0．4115。在任意缺失比例下，非參缺失森林填補(bǔ)算法的優(yōu)勢較明顯。

圖1 不同填補(bǔ)方法在不同缺失比例下的填補(bǔ)效果（NRMSE值）

在不同的缺失比例下，用非參缺失森林填補(bǔ)法對(duì)不同缺失比例下的模擬數(shù)據(jù)集填補(bǔ)10次，均方根誤差的標(biāo)準(zhǔn)差和置信區(qū)間見表1，在5%、10%、20%和30%的缺失比例下，其均方根誤差的標(biāo)準(zhǔn)差分別為0．0006、0．0006、0．0016和0．0312，置信區(qū)間的寬度分別為0．0008、0．008、0．0016和0．0446，說明該算法的穩(wěn)定性強(qiáng)、精確度高。

表1 不同缺失比例下NRMSE均值及標(biāo)準(zhǔn)差（填補(bǔ)10次時(shí)）

2．基于聚類分析的效果評(píng)價(jià)

圖2為四種填補(bǔ)方法在兩個(gè)數(shù)據(jù)集中不同缺失比例下基因功能聚類分析的準(zhǔn)確率（F值）。在不同缺失比例下，不同填補(bǔ)方法對(duì)數(shù)據(jù)集的聚類效果有較大的影響；填補(bǔ)方法上，使用非參缺失森林算法填補(bǔ)數(shù)據(jù)集的聚類效果優(yōu)于其他三種算法。在5%缺失比例的時(shí)候，所有填補(bǔ)方法的F值都高于0．93，聚類效果好；乳腺癌數(shù)據(jù)集在20%缺失比例的時(shí)候，非參缺失森林法、K鄰近填補(bǔ)法、貝葉斯線性回歸法和蒙特卡洛多重填補(bǔ)法的F值依次為0．8819、0．8717、0．7934和0．7501，整體趨勢上和模擬數(shù)據(jù)集中的聚類效果一致。

圖2 不同填補(bǔ)方法在不同缺失比例下基因功能聚類分析的準(zhǔn)確率（F值）

討 論

本文采用不同的方法對(duì)含有缺失值的不同數(shù)據(jù)集進(jìn)行填補(bǔ)，并應(yīng)用于后續(xù)的基因功能聚類分析，通過NRMSE和聚類效果（F值）來評(píng)價(jià)各種填補(bǔ)方法的優(yōu)劣及其適用性，不僅發(fā)展和豐富了基因表達(dá)譜缺失數(shù)據(jù)的填補(bǔ)模型方法，而且為基因表達(dá)譜數(shù)據(jù)分析技術(shù)提供了生物信息學(xué)方法方面的指導(dǎo)。

盡管在不同的數(shù)據(jù)集上依據(jù)不同的指標(biāo)對(duì)各缺失值處理方法進(jìn)行評(píng)價(jià)，結(jié)論會(huì)有細(xì)微的差別，但總體來看，隨著缺失比例的增加，基于非參缺失森林的填補(bǔ)方法優(yōu)勢逐漸顯現(xiàn)，它既提高了缺失估計(jì)的精度和穩(wěn)定性，又可以保留較多的基因信息供后續(xù)的功能分析且具有很高的聚類效能。在缺失比例比較小的情況下，蒙特卡洛多重填補(bǔ)法和K鄰近距離加權(quán)法的填補(bǔ)效果也比較好；從聚類結(jié)果的準(zhǔn)確率來看，如果運(yùn)用不恰當(dāng)?shù)奶钛a(bǔ)方法會(huì)對(duì)后續(xù)表達(dá)譜的研究起誤導(dǎo)性作用，但是直接對(duì)含有缺失數(shù)據(jù)的乳腺癌數(shù)據(jù)進(jìn)行聚類，效果不理想，這也從側(cè)面說明了根據(jù)缺失數(shù)據(jù)集特點(diǎn)選擇正確填補(bǔ)方法的重要性。

本文介紹的非參缺失森林填補(bǔ)方法具有良好的應(yīng)用前景，它對(duì)數(shù)據(jù)集的結(jié)構(gòu)要求較低、實(shí)現(xiàn)簡單，相比于傳統(tǒng)填補(bǔ)算法具有更好的穩(wěn)定性和準(zhǔn)確度，可以保留較多的基因信息供后續(xù)的功能聚類等分析目的。有關(guān)缺失森林程序包的更多擴(kuò)展功能參見missForest程序包說明。本研究結(jié)果是基于較大樣本量且只用于表達(dá)譜數(shù)據(jù)的聚類分析目的，將其推廣到小樣本數(shù)據(jù)和其他分析目的（如差異表達(dá)基因篩選和基因調(diào)控網(wǎng)絡(luò)建立等），可能會(huì)受到一定限制，我們將繼續(xù)進(jìn)行后續(xù)的分析與探討?？傊疚耐ㄟ^不同填補(bǔ)方法的研究，為基因表達(dá)譜數(shù)據(jù)缺失填補(bǔ)策略的建立和缺失填補(bǔ)方法對(duì)基因表達(dá)譜后續(xù)不同分析目的生物學(xué)影響及其程度的評(píng)估，打下了堅(jiān)實(shí)的理論和實(shí)踐基礎(chǔ)。

［1］Nanni L，M ing J，Du Y，et al．M issing value imputation for gene expression data：computational techniques to recovermissing data from available information．American Journal of Medical Genetics，2011，12（5）：498-513．

［2］Troyanskaya O，Cantor M，Sherlock G，etal．M issing value estimation methods for DNA m icroarrays．Bioinformatics，2001，17（6）：520-525．

［3］Kim H，Golub G．M issing value estimation for DNA m icroarray gene expression data：local least squares imputation．Bioinformatics，2005，21（2）：187-198．

［4］Oba S，Sato M，Takemasa I，et al．A Bayesian m issing value estimation method for gene expression profile data．Bioinformatics，2003，volume 19（16）：2088-2096．

［5］J?rnsten R，Wang H，Welsh W，et al．DNA m icroarray data imputation and significance analysis of differential expression．Bioinformatics，2005，21（22）：4155-4161．

［6］Guo X，Alan W，Hong Y．M icroarray m issing data imputation based on a set theoretic framework and biological know ledge．Nucleic Acids Research，2006，34（5）：1608-1619．

［7］Bai F，Liu H．M issing value imputation for m icroarray gene expression data using histone acetylation information．Smart Sensors＆Sensing Technology，2008，9（1）：1-17．

［8］武瑞仙，鄧子兵，譙治蛟，等．利用Monte Carlo技術(shù)模擬研究不同缺失值處理方法對(duì)完全隨機(jī)缺失數(shù)據(jù)的處理效果．中國衛(wèi)生統(tǒng)計(jì)，2015（3）：534-536．

［9］沈琳，胡國清，陳立章，等．缺失森林算法在缺失值填補(bǔ)中的應(yīng)用．中國衛(wèi)生統(tǒng)計(jì)，2014（5）：774-776．［10］Stekhoven D，Bühlmann P．M issForest-non-parametric m issing value imputation for m ixed-type data．Bioinformatics，2012，28（1）：112-118．

［11］Hapfelmeier A，Hothorn T，Riediger C，et al．M ice：multivariate imputation by chained equations in R．International Journal of Biostatistics，2014，45（2）：1-67．

［12］劉熙，王崇駿，葉亮，等．基于最大頻繁項(xiàng)集的層次聚類方法．廣西師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2009，27（3）：105-108．

［13］康茜，李德玉，王素格，等．傳播過程中信號(hào)缺失的層次聚類社區(qū)發(fā)現(xiàn)算法．計(jì)算機(jī)工程與應(yīng)用，2015（9）：201-206．

［14］黃健斌，康劍梅，齊俊杰，等．一種基于同步動(dòng)力學(xué)模型的層次聚類方法．中國科學(xué)（信息科學(xué)），2013（05）：599-610．

［15］楊燕，靳蕃，KAMEL M．聚類有效性評(píng)價(jià)綜述．計(jì)算機(jī)應(yīng)用研究，2008，25（6）：1630-1632．

［16］Gene expression data in estrogen receptor alpha positive breast tumors with and without PIK3CA mutations［http：//www．ncbi．nlm．nih．gov/bioproject/PRJNA128895．

（責(zé)任編輯：劉 壯）

國家自然科學(xué)基金項(xiàng)目（81273178，81573254）

△通信作者：伍亞舟，E-mail：asiawu5＠sina．com