孟繁馳，李書琴，蔡 騁

(西北農(nóng)林科技大學(xué)信息工程學(xué)院，陜西楊凌712100)

0 引言

基因微陣列技術(shù)在生物實(shí)驗(yàn)中已經(jīng)得到了廣泛的應(yīng)用，從微陣列得到的基因表達(dá)值，通常以大規(guī)模矩陣的形式呈現(xiàn)［1］。由于微陣列實(shí)驗(yàn)中各種因素的影響，如雜交失敗、圖像的噪聲、污染等，從微陣列得到的基因矩陣通常含有缺失點(diǎn)［2］。然而，基因下游的實(shí)驗(yàn)，如常用的分類、聚簇方法等，通常需要完整的矩陣作為輸入。鑒于此，不同類型的基因微陣列缺失點(diǎn)的重建方法陸續(xù)產(chǎn)生。目前常用的基因微陣列缺失點(diǎn)的重建方法包括奇異值分解SVD［3］、K最近鄰KNN［3］、貝葉斯主成分分析BPCA［4］、最小二乘法LSimpute［5］、局部最小二乘法 LLSimpute［6］等。在以上各種方法中，LLSimpute被實(shí)驗(yàn)證明可以在各種類型的基因數(shù)據(jù)集上取得較高的重建精確度［1，6-8］。

1 LLSimpute方法

LLSimpute是Hyunsoo Kim等人在2005年提出的方法［6］。該方法利用基因與基因之間的線性關(guān)系來重建目標(biāo)基因中的缺失點(diǎn)。以一個(gè)m×n的基因微陣列矩陣為例，簡要介紹一下LLSimpute的流程:

(1)從其他m-1個(gè)基因行中找出與目標(biāo)基因 (即當(dāng)前需要重建的含缺失點(diǎn)的一行基因)g1距離最近的k個(gè)基因，組成列向量 (g1，gs1，…，gsk)T。

(2)找出目標(biāo)基因中的缺失點(diǎn)g1×q和非缺失點(diǎn)g1×(n-q)。其中q表示該行中有q個(gè)缺失點(diǎn)。

(3)將缺失的點(diǎn)g1×q組成的向量表示為α，將非缺失的點(diǎn)g1×(n-q)組成的向量表示為wT。從gs1到gsk中找出與g1中的缺失點(diǎn)對應(yīng)的列上的基因，組成矩陣 (或向量)Bk×q，從gs1到gsk中找出與g1中的非缺失點(diǎn)對應(yīng)的列上的基因，組成矩陣Ak×(n-q)。假設(shè)目標(biāo)基因上有兩個(gè)缺失點(diǎn)，則形成的矩陣如式 (1)所示

(4)在以上矩陣中，wT、B、A是已知的，α表示目標(biāo)基因中待重建的缺失點(diǎn)組成的向量。由wT和A之間的線性關(guān)系，求系數(shù)向量x，可以看作一個(gè)最小二乘法問題 (2)

其中(AT)是AT的偽逆陣 (Pseudoinverse)。從而

實(shí)驗(yàn)表明，在k選擇合適的時(shí)候，LLSimpute可以獲得非常好的重建效果。［6］中，為了找出最優(yōu)的k值，作者采取了如下的方法:首先，用每行的均值將該行中的缺失點(diǎn)填充，得到一個(gè)完整的矩陣。然后在該完整的矩陣當(dāng)中人為地去掉一些點(diǎn)，構(gòu)成一個(gè)人造的含缺失點(diǎn)的矩陣。用不同數(shù)量的k個(gè)鄰居，使用LLSimpute去重建，得到重建出的缺失點(diǎn)的數(shù)值。由于這些缺失點(diǎn)是人為制造的，所以它們的真實(shí)值是已知的，這樣，不同k對應(yīng)的重建錯誤率是可以計(jì)算的。最后，取對應(yīng)重建錯誤率最低的k，在基因微陣列矩陣中作為選取的鄰居數(shù)。

2 矩陣填充

矩陣填充 (matrix completion，MC)是 Candès等人在2009年提出的方法 ［9］。在 ［9］中，作者證明了一個(gè)含有缺失點(diǎn)的矩陣，在低秩的情況下，如果其中所含的非缺失點(diǎn)足夠多，則完整的矩陣是可以在一定的概率下恢復(fù)出來的。因此，矩陣填充可以認(rèn)為是一個(gè) (4)所示的優(yōu)化問題

其中M是觀測到的，含有缺失點(diǎn)的矩陣，Ω是非缺失點(diǎn)的元素對應(yīng)的下標(biāo)組成的集合，X是待重建的完整矩陣。但 (4)是一個(gè)NP難問題。然而矩陣X的秩與矩陣X的奇異值中非零元素的個(gè)數(shù)是一樣的，所以［9］中，X的秩被X的核范數(shù) (nuclear norm)所代替，上式變?yōu)?(5)所示的優(yōu)化問題

其中‖·‖*表示核范數(shù)，等于所有奇異值的和。式(4)和式 (5)的區(qū)別在于，式 (4)中X的秩對應(yīng)的是其奇異值組成的向量的0范數(shù)，而式 (5)中X的核范數(shù)對應(yīng)其奇異值組成的向量的1范數(shù)。式 (5)所對應(yīng)的優(yōu)化問題是凸的，具有唯一的解。

矩陣填充已經(jīng)被成功用在Netflix視頻推薦系統(tǒng)［10］、視頻去噪上［11］。林宙辰等人提出的非精確增廣拉格朗日乘子(inexact augmented lagrange multiplier，IALM)［12］是解決式(5)所示的優(yōu)化問題的一個(gè)有效方法，在重建精確度和運(yùn)行速度上均占有很大優(yōu)勢。

IALM采用拉格朗日乘子解決 (5)所示的核范數(shù)凸優(yōu)化問題，這里，矩陣填充問題被描述為式 (6)所示的一個(gè)優(yōu)化問題

式中:A——待重建的矩陣，D——觀測到的有缺失點(diǎn)的矩陣，πΩ——一個(gè)的線性變換。IALM在每次迭代中求奇異值分解，每次取較大的幾個(gè)奇異值來重建矩陣A和E，從而使A趨向于低秩。詳細(xì)的迭代過程，可參考［12］。本文中，將使用IALM進(jìn)行矩陣填充，得到完整基因微陣列矩陣的方法稱為MC方法。該方法和其他方法的對比實(shí)驗(yàn)結(jié)果將在第4章中給出。

3 矩陣填充對LLSimpute的改進(jìn)

由LLSimpute中的第1步可以看出來，與目標(biāo)基因?qū)?yīng)的k最近鄰基因行必須是完整的，才可以構(gòu)成式中的A矩陣和B矩陣。同樣，在自適應(yīng)地尋找最優(yōu)的k時(shí)，也需要完整的矩陣。然而實(shí)際的基因微陣列矩陣往往是多行都含有缺失點(diǎn)。一種解決辦法是在找k最近鄰的時(shí)候，忽略掉含有缺失點(diǎn)的基因，只在完整的基因中尋找鄰居，這種辦法在缺失率比較低的時(shí)候是可行的，但是當(dāng)缺失率較高時(shí)，大多數(shù)基因都是不完整的，按照這個(gè)方法，會帶來鄰居數(shù)量過少的問題。實(shí)驗(yàn)也表明，LLSimpute對k選取是非常敏感的，若只在完整的行中尋找鄰居，在缺失率大于等于5%的時(shí)候，LLSimpute就會因可選取的鄰居數(shù)量過少，而使重建錯誤率大大升高。

由此可以看出，LLSimpute需要一個(gè)事先填充好的“中間矩陣”才可以構(gòu)造A矩陣和B矩陣，以及尋找最優(yōu)的k。對于這個(gè)“中間矩陣”，文獻(xiàn)［6］中采用行均值的方法得到，即在含缺失點(diǎn)的矩陣當(dāng)中，用每行的均值來代替缺失點(diǎn)，得到完整的矩陣。在本文的第2章中，已經(jīng)提到了用矩陣填充得到完整矩陣的方法。矩陣填充得到的結(jié)果，可以作為最終重建后的微陣列矩陣，也可以作為一個(gè)“中間矩陣”，繼續(xù)使用LLSimpute來進(jìn)行重建。這樣，即將LLSimpute中的行均值填充的過程，替換為矩陣填充的過程。在本文中，將這種改進(jìn)的LLSimpute稱為MC_LLS方法。

用行均值填充可以構(gòu)成一個(gè)完整的矩陣，但行均值只利用了每個(gè)目標(biāo)基因本行的信息。MC_LLS方法中的MC過程是以矩陣的低秩為優(yōu)化目標(biāo)，利用的是矩陣的全局特性，而LLSimpute本身，是基于矩陣的局部特征的，所以MC_LLS方法，實(shí)際上同時(shí)利用了矩陣的全局特征和局部特征。MC_LLS方法的實(shí)驗(yàn)結(jié)果，也將在第4章中給出。

4 實(shí)驗(yàn)與結(jié)果

對目前常用的幾種方法 (KNN、BPCA、LLSimpute)，和本文提出的MC、MC_LLS方法進(jìn)行對比試驗(yàn)，并給出結(jié)果的分析。

4.1 數(shù)據(jù)集

試驗(yàn)在4個(gè)公開的真實(shí)基因微陣列數(shù)據(jù)集上進(jìn)行。前兩個(gè)數(shù)據(jù)都屬于時(shí)間序列。其中第一個(gè)數(shù)據(jù)集來自［13］中的酵母細(xì)胞的CDC15序列和CDC28序列，將其稱為CDC15_28。第二個(gè)來自 ［13］中的酵母細(xì)胞的alpha序列，稱為SP_APLHA。第三個(gè)數(shù)據(jù)集來自文獻(xiàn)［14］中的人類癌細(xì)胞株，是一個(gè)非時(shí)間序列，稱為NCI60。第四個(gè)數(shù)據(jù)集來自［15］中的淋巴瘤細(xì)胞，也是一個(gè)非時(shí)間序列，稱為lymphoma。其中CDC15_28在綜述［1］中也被用于對各種方法進(jìn)行比較，SP_APLHA在 ［6］中被用于檢驗(yàn)LLSimpute的重建精確度，NCI60與綜述［1］中的NTS來自相同的細(xì)胞株，lymphoma在LSimpute［5］中也被采用。各個(gè)數(shù)據(jù)集的總行列數(shù)、完整的行列數(shù) (即除去含缺失點(diǎn)的行之后的大小)以及類型可見表1。

表1 實(shí)驗(yàn)采用的數(shù)據(jù)集

這些數(shù)據(jù)集原本都是不完整的，為了評價(jià)重建后的錯誤率，需要完整的微陣列矩陣作為參考。為此，將矩陣中不完整的行去掉，只保留完整的行 (如表1中的第三列)。再在這些完整的行中按不同的比例隨機(jī)去除一些已知點(diǎn)，作為缺失點(diǎn)。同樣的方法也在 ［1，6］中被采用。這樣，可以同時(shí)得到完整的微陣列矩陣和它們對應(yīng)的不同缺失率下的缺失矩陣。對不同缺失率的矩陣進(jìn)行重建，將重建后的矩陣和完整的矩陣對比，可以得出重建誤差。

4.2 參數(shù)選擇

對于KNN方法，在CDC15_28和SP_ALPHA數(shù)據(jù)集上，鄰居數(shù)k值設(shè)為10;在NCI60和lymphoma數(shù)據(jù)集上，k設(shè)為5，這兩個(gè)k值處于KNN最優(yōu)的參數(shù)范圍區(qū)間內(nèi)［1］。BPCA的參數(shù)設(shè)為其默認(rèn)值，即微陣列矩陣的列數(shù)減去1。LLSimpute中的參數(shù)k使用 ［6］中模擬缺失點(diǎn)篩選k的方法得到的值。MC_LLS中的k也采用模擬缺失點(diǎn)的方法得到。需要說明的是，對于LLSimpute和MC_LLS，在尋找鄰居時(shí)，若缺失率足夠低 (此時(shí)完整的行數(shù)相對比較多)，則只從完整的基因行中尋找，而不是從用行均值或MC填充好的“中間矩陣”中尋找，以有效利用矩陣本身的真實(shí)值。

4.3 重建評價(jià)標(biāo)準(zhǔn)

目前對于基因重建效果的評價(jià)，多采用的是標(biāo)準(zhǔn)化根均方差 (normalized root mean squared error)［1，4，6］，稱為NRMSE，定義如下

式中:yj——基因第j位上原始的值，j——該位置上重建后的值，N——要重建的基因點(diǎn)的總數(shù)，即缺失點(diǎn)的個(gè)數(shù)，σy——N個(gè)原始基因點(diǎn)的標(biāo)準(zhǔn)差 (standard deviation)。NRMSE值越低，表明重建錯誤率越低。

4.4 重建結(jié)果

對4.1中的4個(gè)數(shù)據(jù)集，按照從1%到25%的之間的8種不同的缺失率，構(gòu)造缺失點(diǎn)。對于每種數(shù)據(jù)集，不同缺失率下的重建都進(jìn)行10次重復(fù)實(shí)驗(yàn)，取10次標(biāo)準(zhǔn)化根均方差的平均值作為最終的重建錯誤率。圖1—圖4分別給出了不同方法在CDC15_28、SP_ALPHA、NCI60和lymphoma數(shù)據(jù)集上的結(jié)果。

圖1 在數(shù)據(jù)集CDC15_28上的重建錯誤率

圖1中，KNN的曲線在較高的缺失率下被截?cái)?，表示在這些錯誤率下，基因的每行都存在缺失點(diǎn)，從而無法從完整的行中找到最近鄰。在缺失率小于10%的情況下，MC的重建錯誤率明顯低于KNN、BPCA和LLSimpute。但是隨著缺失率的升高，MC的錯誤率也開始變大。然而MC_LLS表現(xiàn)很穩(wěn)定，在任何缺失率下，重建錯誤率都是最低的。

圖2 在數(shù)據(jù)集SP_ALPHA上的重建錯誤率

由圖2可以看出，在數(shù)據(jù)集SP_ALPHA下，KNN的表現(xiàn)仍然是最差的，重建錯誤率始終高于0.9。LLSimpute在缺失率較低的情況下有較低的錯誤率，但當(dāng)缺失率大于5%時(shí)，錯誤率明顯提高。MC在缺失率為8%的時(shí)候，錯誤率小于LLSimpute，在其他缺失率下均高于LLSimpute、BPCA和MC_LLS。但是總體上，仍然是MC_LLS取得了最低的重建錯誤率。

圖3 在數(shù)據(jù)集NCI60上的重建錯誤率

NCI60數(shù)據(jù)集是一個(gè)非時(shí)間序列，具有很強(qiáng)的局部相關(guān)性，所以圖3中基于數(shù)據(jù)全局特征的BPCA和MC在這里的重建錯誤率明顯比LLS、MC_LLS高出很多。而基于局部相關(guān)性的KNN，在這里的表現(xiàn)雖然仍然比較差，但是已經(jīng)取得了比時(shí)間序列下更低的錯誤率。對于MC_LLS，在大多數(shù)的缺失率下，仍然具有最低的重建錯誤率。

圖4 在數(shù)據(jù)集lymphoma上的重建錯誤率

與數(shù)據(jù)集NCI60類似，圖4中的數(shù)據(jù)集lymphoma也具有較強(qiáng)的局部相關(guān)性。LLSimpute和BPCA的表現(xiàn)非常接近，KNN的結(jié)果仍然是最差的。雖然MC本身的NRMSE高出LLSimpute和BPCA，但MC_LLS仍然在除1%之外的其他所有缺失率下取得了最低的重建錯誤率。

由以上實(shí)驗(yàn)結(jié)果可以得出結(jié)論，基于矩陣全局特征的MC方法在局部相關(guān)性較弱的數(shù)據(jù)集上具有優(yōu)勢，但是不適用于局部相關(guān)性很強(qiáng)的數(shù)據(jù)集。而MC_LLS由于結(jié)合了矩陣的整體相關(guān)性和局部相關(guān)性，所以在任何類型的數(shù)據(jù)集上都表現(xiàn)出了最好的重建結(jié)果。

表2所示的是各種方法在CDC15_28數(shù)據(jù)集上，不同缺失率下的運(yùn)算時(shí)間，單位是秒。每個(gè)時(shí)間都是在進(jìn)行了10次重復(fù)實(shí)驗(yàn)后求得的平均值。實(shí)驗(yàn)在一臺裝有4核Intel Xeon E5504 2.0GHz處理器的計(jì)算機(jī)上進(jìn)行，采用Matlab R2011a。

表2 各種方法在數(shù)據(jù)集CDC15_28上的運(yùn)算時(shí)間 (秒)

雖然矩陣填充中的奇異值分解通常是一項(xiàng)比較耗時(shí)的計(jì)算，但由于IALM使用了部分奇異值分解［12］，所需時(shí)間大大減少，所以MC方法在運(yùn)算時(shí)間上具有明顯的優(yōu)勢。而相對于LLSimpute，MC_LLS在時(shí)間開銷增加不明顯的情況下，卻明顯降低了重建錯誤率。

5 結(jié)束語

使用非精確增廣拉格朗日乘子 (IALM)解決了矩陣填充 (MC)問題，實(shí)現(xiàn)了基因微陣列矩陣的核范數(shù)凸優(yōu)化。含有缺失點(diǎn)的微陣列矩陣，在進(jìn)行矩陣填充后，可以作為最終的重建結(jié)果，并且矩陣填充得到的結(jié)果作為中間矩陣，還可以替換LLSimpute中的行均值矩陣。在四個(gè)常用的基因微陣列數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果表明矩陣填充方法 (MC)在某些時(shí)間序列微陣列矩陣的部分缺失率下可以取得比KNN、LLSimpute和BPCA更高的重建精確度，并且計(jì)算時(shí)間具有明顯優(yōu)勢。矩陣填充和LLSimpute結(jié)合的方法 (MC_LLS)，在實(shí)驗(yàn)的四個(gè)數(shù)據(jù)集中的幾乎所有缺失率下，取得了最高的重建精確度，并且相對LLSimpute時(shí)間增加不明顯。

［1］Bras L P，Menezes J C.Dealing with gene expression missing data［J］.Syst Biol(Stevenage)，2006，153(3):105-19.

［2］Liew Alan Wee-Chung，Law Ngai-Fong，Yan Hong.Missing value imputation for gene expression data:Computational techniques to recover missing data from available information［J］.Briefings in Bioinformatics，2011，12(5):498-513.

［3］Troyanskaya Olga，Cantor Michael，Sherlock Gavin，et al.Missing value estimation methods for DNA microarrays［J］.Bioinformatics，2001，17(6):520.

［4］Oba Shigeyuki，Sato Masa-aki，Takemasa Ichiro，et al.A bayesian missing value estimation method for gene expression profile data［J］.Bioinformatics，2003，19(16):2088.

［5］B Trond Hellen，Dysvik Bjarte，Jonassen Inge.LSimpute:Accurate estimation of missing values in microarray data with least squares methods［J］.Nucleic Acids Research，2004，32(3):e34.

［6］Kim Hyunsoo，Golub Gene H，Park Haesun.Missing value estimation for DNA microarray gene expression data:Local least squares imputation［J］.Bioinformatics，2005，21(2):187.

［7］Brock Guy N，Shaffer John R，Blakesley Richard E，et al.Which missing value imputation method to use in expression profiles:A comparative study and two selection schemes［J］.BMCbioinformatics，2008，9(1):12.

［8］Magalie C，Alain M，Ga lle L.Comparative analysis of missing value imputation methods to improve clustering and interpretation of microarray experiments［J］.BMCgenomics，2010，11(1):15.

［9］Candès Emmanuel J，Recht Benjamin.Exact matrix completion via convex optimization［J］.Foundations of Computational Mathematics，2009，9(6):717-772.

［10］Bennett James，Lanning Stan.The netflix prize［R］.California，USA:KDDCup，2007.

［11］JI Hui，LIU Chaoqiang，SHEN Zuowei，et al.Robust video denoising using low rank matrix completion［C］//IEEE Conference on Computer Vision and Pattern Recognition，2010.

［12］LIN Zhouchen，CHEN Minming，MA Yi.The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices［J］.Arxiv preprint arXiv:1009.5055，2010.

［13］Spellman Paul T，Sherlock Gavin，ZHANG Michael Q，et al.Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization［J］.Molecular biology of the cell，1998，9(12):3273.

［14］Scherf Uwe，Ross Douglas T，Waltham Mark，et al.A gene expression database for the molecular pharmacology of cancer［J］.Nature genetics，2000，24(3):236-244.

［15］Alizadeh A Alizadeh，Eisen Michael B，Davis R.Eric，et al.Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling［J］.Nature，2000，403(6769):503-511.