基因表達譜的非參缺失森林填補算法研究*

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

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

·方法介紹·

基因表達譜的非參缺失森林填補算法研究*

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

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

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

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

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

理論與算法

1．非參缺失森林的填補方法

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

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

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

（1）首先用應變量Yobs和自變量Xobs擬合一個隨機森林；

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

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

2．缺失填補的效果評價

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

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

3．聚類的效能評價

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

（1）初始時共有N類，每個類有一個對象構(gòu)成。令序號m＝0，l（m）＝0。

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

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

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

（5）重復（2）～（4）步，直到所有對象合并成一個類為止。

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

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

數(shù)據(jù)集

1．數(shù)據(jù)來源

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

2．統(tǒng)計分析

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

結(jié) 果

1．基于均方根的填補效果評價

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

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

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

表1 不同缺失比例下NRMSE均值及標準差（填補10次時）

2．基于聚類分析的效果評價

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

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

討 論

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

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

本文介紹的非參缺失森林填補方法具有良好的應用前景，它對數(shù)據(jù)集的結(jié)構(gòu)要求較低、實現(xiàn)簡單，相比于傳統(tǒng)填補算法具有更好的穩(wěn)定性和準確度，可以保留較多的基因信息供后續(xù)的功能聚類等分析目的。有關(guān)缺失森林程序包的更多擴展功能參見missForest程序包說明。本研究結(jié)果是基于較大樣本量且只用于表達譜數(shù)據(jù)的聚類分析目的，將其推廣到小樣本數(shù)據(jù)和其他分析目的（如差異表達基因篩選和基因調(diào)控網(wǎng)絡建立等），可能會受到一定限制，我們將繼續(xù)進行后續(xù)的分析與探討?？傊?，本文通過不同填補方法的研究，為基因表達譜數(shù)據(jù)缺失填補策略的建立和缺失填補方法對基因表達譜后續(xù)不同分析目的生物學影響及其程度的評估，打下了堅實的理論和實踐基礎(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ù)模擬研究不同缺失值處理方法對完全隨機缺失數(shù)據(jù)的處理效果．中國衛(wèi)生統(tǒng)計，2015（3）：534-536．

［9］沈琳，胡國清，陳立章，等．缺失森林算法在缺失值填補中的應用．中國衛(wèi)生統(tǒng)計，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］劉熙，王崇駿，葉亮，等．基于最大頻繁項集的層次聚類方法．廣西師范大學學報（自然科學版），2009，27（3）：105-108．

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

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

［15］楊燕，靳蕃，KAMEL M．聚類有效性評價綜述．計算機應用研究，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．

（責任編輯：劉 壯）

國家自然科學基金項目（81273178，81573254）

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