基于模擬退火優(yōu)化雙聚類的基因數(shù)據(jù)填補(bǔ)方法

朱 嫻 楊 明 馬 衛(wèi) 朱 俊

1(南京理工大學(xué)紫金學(xué)院計(jì)算機(jī)學(xué)院 江蘇 南京 210046) 2(南京旅游職業(yè)學(xué)院酒店管理學(xué)院 江蘇 南京 211100)

基于模擬退火優(yōu)化雙聚類的基因數(shù)據(jù)填補(bǔ)方法

朱 嫻1楊 明1馬 衛(wèi)2朱 俊1

1(南京理工大學(xué)紫金學(xué)院計(jì)算機(jī)學(xué)院 江蘇 南京 210046)2(南京旅游職業(yè)學(xué)院酒店管理學(xué)院 江蘇 南京 211100)

基因表達(dá)數(shù)據(jù)是由DNA微陣列實(shí)驗(yàn)產(chǎn)生的大規(guī)模矩陣，能有效地提取生物學(xué)信息，由于受到實(shí)驗(yàn)條件限制，基因表達(dá)數(shù)據(jù)往往存在缺失值，需要進(jìn)行缺失數(shù)據(jù)的填補(bǔ)。傳統(tǒng)的缺失數(shù)據(jù)填補(bǔ)方法是基于基因表達(dá)數(shù)據(jù)的單一特征，未充分考慮數(shù)據(jù)矩陣間的相關(guān)性。針對(duì)雙聚類均方殘值越小基因表達(dá)數(shù)據(jù)相關(guān)性越高這一特性進(jìn)行研究，提出一種基于模擬退火優(yōu)化雙聚類的缺失數(shù)據(jù)填補(bǔ)方法(bi-SA)，采用模擬退火法確定最優(yōu)雙聚類，從而實(shí)現(xiàn)缺失數(shù)據(jù)的最有效填補(bǔ)。四組真實(shí)基因表達(dá)數(shù)據(jù)實(shí)驗(yàn)表明，bi-SA方法能夠獲得較高的填補(bǔ)準(zhǔn)確性。

基因表達(dá)數(shù)據(jù) 缺失數(shù)據(jù)填補(bǔ) 模擬退火法 雙聚類

0 引 言

基因芯片技術(shù)[1]，又稱為DNA微陣列技術(shù)，是一種高通量的分子生物學(xué)手段，能夠?qū)Υ罅康幕蜻M(jìn)行快速測(cè)量，同時(shí)獲得成千上萬(wàn)個(gè)基因在不同條件下的表達(dá)水平?；蛐酒ぷ髟韀2]是對(duì)兩個(gè)互補(bǔ)DNA分子雜交信號(hào)進(jìn)行檢測(cè)分析，獲取待測(cè)生物樣本的基因表達(dá)豐度，在生物學(xué)和醫(yī)學(xué)研究等方面獲得了廣泛的應(yīng)用，為探索生命科學(xué)提供了強(qiáng)有力的工具。從基因芯片實(shí)驗(yàn)產(chǎn)生的基因表達(dá)數(shù)據(jù)通常是一個(gè)大形式的矩陣，不可避免地受到技術(shù)、實(shí)驗(yàn)條件的限制而產(chǎn)生數(shù)據(jù)缺失問(wèn)題?；虮磉_(dá)數(shù)據(jù)分析包括聚類、分類、特征值提取等方法，一般分析需要輸入完整的數(shù)據(jù)矩陣，對(duì)于具有缺失值的不完備基因表達(dá)數(shù)據(jù)矩陣不能直接處理。最初的解決方法是重復(fù)實(shí)驗(yàn)找出缺失點(diǎn)或者刪除包含缺失值的基因，但是這兩種方法實(shí)用價(jià)值不高。隨后采用0值或行均值[3]填補(bǔ)缺失值的方法雖然簡(jiǎn)單，但是處理數(shù)據(jù)的方式過(guò)于粗糙、誤差率高，對(duì)后續(xù)實(shí)驗(yàn)結(jié)果產(chǎn)生較大影響，在底層的基因表達(dá)數(shù)據(jù)挖掘分析中被證明是缺少有用價(jià)值的。

數(shù)據(jù)缺失問(wèn)題在各種應(yīng)用研究領(lǐng)域都會(huì)出現(xiàn)，相關(guān)的缺失數(shù)據(jù)的填補(bǔ)方法研究歷史比較悠久，針對(duì)簡(jiǎn)單填補(bǔ)法的不足，提出多重填補(bǔ)法MI(multiple imputation)[4-5]。通過(guò)對(duì)數(shù)據(jù)的填補(bǔ)、分析、合并，將數(shù)據(jù)進(jìn)行綜合分析從而統(tǒng)計(jì)推斷出數(shù)據(jù)。近幾年，基因表達(dá)數(shù)據(jù)缺失值填補(bǔ)方法研究陸續(xù)提出，主要集中在混合全局的局部填補(bǔ)方法和雙聚類填補(bǔ)方法。

混合全局的局部填補(bǔ)方法通過(guò)全局和局部多種方法的組合估計(jì)缺失數(shù)據(jù)，利用統(tǒng)計(jì)分析的方法捕獲數(shù)據(jù)中全局和局部相關(guān)性。既能克服全局填補(bǔ)算法產(chǎn)生過(guò)擬合問(wèn)題，又能彌補(bǔ)局部填補(bǔ)方法中基因表達(dá)數(shù)據(jù)表達(dá)豐度分布不均勻的缺點(diǎn)。Jornsten等[6]通過(guò)加權(quán)和獲得基因表達(dá)數(shù)據(jù)的估計(jì)值，Pan等[7]在此基礎(chǔ)上對(duì)權(quán)值求解進(jìn)行改進(jìn)。Wang等[8]通過(guò)皮爾遜系數(shù)獲得目標(biāo)基因，通過(guò)局部最小二乘法(LLS)和收縮估計(jì)法獲得缺失值。局部最小二乘法(LLS)是最常用的局部基因表示矩陣缺失點(diǎn)估計(jì)方法，采用多元回歸模型，具有計(jì)算復(fù)雜度低、估計(jì)精度高等優(yōu)點(diǎn)。Li等[9]和He等[10]均采用LLS法，結(jié)合全局特征的貝葉斯主元分析法，利用主軸向量的線性組合估計(jì)基因表達(dá)數(shù)據(jù)中的缺失值。Cheng等[11]將雙聚類應(yīng)用在基因表達(dá)數(shù)據(jù)矩陣中，目的是找出具有相近表達(dá)水平的子矩陣。Ji等[12]將雙聚類和數(shù)據(jù)矩陣缺失值的估計(jì)結(jié)合，通過(guò)雙聚類獲得局部結(jié)構(gòu)，充分利用子矩陣與缺失點(diǎn)相關(guān)性最強(qiáng)的信息，降低估算錯(cuò)誤率。bi-BPCA[13]方法將k近鄰算法(KNN)、貝葉斯主元分析法(PCA)與雙聚類相結(jié)合，克服了貝葉斯主元分析法不能處理數(shù)據(jù)局部結(jié)構(gòu)的缺陷。KNN填補(bǔ)算法[14]是較為流行的一種局部方法，通過(guò)計(jì)算目標(biāo)基因與其他具有完全值基因間的距離選擇最近鄰居基因，能較好地描述基因之間的相關(guān)程度。BIA[15]方法利用雙聚類均方殘值越小數(shù)據(jù)相關(guān)性越高的特點(diǎn)來(lái)預(yù)測(cè)基因表達(dá)數(shù)據(jù)缺失值，并進(jìn)行相關(guān)定理證明雙聚類的均方殘值越低，預(yù)測(cè)結(jié)果的準(zhǔn)確性越高，但BIA方法是在已知雙聚類基礎(chǔ)上產(chǎn)生，對(duì)雙聚類的產(chǎn)生未做進(jìn)一步的研究。本文提出的一種基于模擬退火優(yōu)化雙聚類的缺失數(shù)據(jù)填補(bǔ)方法(bi-SA)，同樣是基于雙聚類原理，能較好地處理基因表達(dá)數(shù)據(jù)的相關(guān)性有效地填補(bǔ)缺失值。

1 雙聚類基本原理

令基因表達(dá)數(shù)據(jù)矩陣為A，aij為表達(dá)矩陣A的元素值，表示在第j列條件下第i行基因的表達(dá)值。在理想狀態(tài)下，雙聚類的元素值aij可以定義為行均值加上列均值減去矩陣均值。

aij=aiJ+aIj-aIJ

(1)

式中:I為基因集合的子集，J為條件集合的子集。然而基因數(shù)據(jù)中存在噪聲，雙聚類并不一定是理想的，因此定義殘差量化一個(gè)元素的真實(shí)值與相關(guān)的行均值、列均值、矩陣均值的期望。非理想狀態(tài)下元素值aij的殘差定義為：

RSIJ(i,j)=aij-aiJ-aIj+aIJ

(2)

為了評(píng)價(jià)雙聚類在噪聲影響下的質(zhì)量，Cheng等[11]在殘差的基礎(chǔ)上定義了均方殘值H(I,J)為雙聚類行與列的相關(guān)性打分。均方殘值表示為：

(3)

如果滿足H(I,J)≤δ，且δ≥0，則子矩陣AIJ叫作δ-雙聚類。

(4)

令x為子矩陣AIJ中的一個(gè)缺失值，行數(shù)為m，列數(shù)為n。其中p為缺失值x所在的行號(hào)，q為缺失元素x所在的列號(hào)，G={1,2,…,p-1,p+1,…,m}C={1,2,…,q-1,q+1,…,n}為非缺失值集合，則缺失值x的最優(yōu)值[15]為：

(5)

2 基于模擬退火的雙聚類缺失數(shù)據(jù)方法

2.1 k近鄰法

k近鄰法是基因表達(dá)數(shù)據(jù)缺失值填補(bǔ)最常用方法之一，它是尋找與目標(biāo)基因距離最近的k個(gè)具有完全值的基因。本文采用歐氏距離計(jì)算基因之間的距離，表示為：

(6)

式中:at表示包含缺失值的目標(biāo)向量，as表示目標(biāo)向量外的任意完全向量。雖然k近鄰算法簡(jiǎn)單使用廣泛，但是基因表達(dá)數(shù)據(jù)相關(guān)性不能有效利用。為了更好地利用基因表達(dá)數(shù)據(jù)的特征，本文采用k近鄰法在行和列兩個(gè)方向上同時(shí)獲取，得到體積較小的矩陣作為雙聚類初始種子S(I,J)。

2.2 模擬退火法

模擬退火法SA(simulated annealing)[16]基本思想是從某一較高溫度開始，伴隨著溫度下降在解空間隨機(jī)尋找目標(biāo)函數(shù)的全局最優(yōu)解，從而有效避免陷入局部最優(yōu)的優(yōu)化算法。模擬退火法常用于組合優(yōu)化，在雙聚類優(yōu)化過(guò)程中，增加雙聚類初始種子S(I,J)的行和列，產(chǎn)生新雙聚類種子S′(I,J)，采用式(3)計(jì)算初始雙聚類的均方殘值，獲得均方殘值增量公式為：

ΔE=H(S′)-H(S)

(7)

當(dāng)ΔE<0，則接受種子S′作為當(dāng)前新的雙聚類，否則計(jì)算狀態(tài)接受函數(shù)。

f=exp(-ΔE/T)

(8)

當(dāng)概率值f>rand(0,1)，則接受種子S′(I,J)作為當(dāng)前新的雙聚類。

2.3 bi-SA方法描述

bi-SA步驟如下所述。

1) 找出基因表達(dá)數(shù)據(jù)集中缺失值x，計(jì)算x所在行的行均值，作為進(jìn)行x初始值。

2) 利用式(6)產(chǎn)生雙聚類初始種子。

3) 模擬退火法優(yōu)化雙聚類。

a) 初始化初始溫度、降溫比率、終止溫度；

b) 利用式(3)計(jì)算雙聚類種子的均方殘值；

c) 增加雙聚類種子的行和列，產(chǎn)生新的雙聚類種子，并計(jì)算均方殘值和均方殘值增量ΔE；

d) 判斷增量ΔE，若增量<0，則接受新解為當(dāng)前解，否則以概率exp(-ΔE/T)接受作為新的當(dāng)前解；

e) 如果滿足終止條件，輸出當(dāng)前缺失值x所在的最優(yōu)雙聚類；

f) 逐漸減少溫度，然后轉(zhuǎn)b步。

4) 利用式(5)計(jì)算最優(yōu)雙聚類缺失值x。

從上述步驟可知bi-SA方法采用k近鄰法產(chǎn)生含有缺失值x的雙聚類種子，根據(jù)模擬退火法降溫比率0.9優(yōu)化目標(biāo)雙聚類，計(jì)算缺失值x。bi-SA方法充分利用基因表達(dá)數(shù)據(jù)的局部特性和雙聚類的生物特性，產(chǎn)生最優(yōu)雙聚類時(shí)保證了缺失值的填補(bǔ)質(zhì)量，很好地反映出真實(shí)缺失數(shù)據(jù)信息。

3 實(shí)驗(yàn)結(jié)果與分析

本文采用四組數(shù)據(jù)Infection[17]、Ronen[18]、Ogawa[19]、Yoshi[20]，如表1所示。

表1 四種實(shí)驗(yàn)數(shù)據(jù)集

其中，Infection和Ronen為時(shí)間序列數(shù)據(jù)集，Ogawa為非時(shí)間序列數(shù)據(jù)集，Yoshi為混合序列數(shù)據(jù)集。為了驗(yàn)證bi-SA填補(bǔ)算法，與雙聚類貝葉斯主成分分析法(bi-BPCA)[13]、貝葉斯主成分分析法(BPCA)、局部最小二乘法(LLS)和雙聚類迭代最小二乘法法(Li-iLs)進(jìn)行對(duì)比實(shí)驗(yàn)。算法編程實(shí)現(xiàn)采用Matlab 2012b。實(shí)驗(yàn)采用在完整行列中按不同比例模擬缺失點(diǎn)的方法，在完整的行列中按一定的比例隨機(jī)去掉數(shù)據(jù)作為缺失值，評(píng)價(jià)采用標(biāo)準(zhǔn)化均方根誤差，表示為：

(9)

圖1為bi-SA方法與四種方法在四組數(shù)據(jù)中的比較。其中：橫軸為數(shù)據(jù)的缺失率，縱軸為標(biāo)準(zhǔn)化均方根誤差。(a)采用Infection數(shù)據(jù)集，當(dāng)缺失率為1%到20%時(shí)，bi-SA方法的標(biāo)準(zhǔn)化均方根誤差最低，但是當(dāng)缺失率高于25%時(shí)，bi-SA方法的錯(cuò)誤率較高，當(dāng)缺失率為30%時(shí)波動(dòng)最大。(b)采用Ronen數(shù)據(jù)集，缺失率為25%時(shí)有較大的波動(dòng)，其他缺失率的錯(cuò)誤率均優(yōu)于其他四種方法。(c)采用Ogawa數(shù)據(jù)集，與Ronen數(shù)據(jù)集類似，僅當(dāng)缺失率為25%時(shí)，錯(cuò)誤率略高于四種算法但無(wú)較大波動(dòng)。(d)采用Yoshi數(shù)據(jù)集，當(dāng)缺失率為20%時(shí)，錯(cuò)誤率小于其他四種方法，其他的缺失率均較高且具有較大的波動(dòng)。

圖1 不同數(shù)據(jù)集比較

通過(guò)四組數(shù)據(jù)比較可以看出bi-SA方法在缺失率為20%時(shí)的錯(cuò)誤率最低，當(dāng)大于20%時(shí)錯(cuò)誤率均明顯的升高。bi-SA方法采用模擬退火法對(duì)雙聚類進(jìn)行搜索，得到的結(jié)果具有概率突跳性，錯(cuò)誤率的高低與雙聚類的容量和均方殘值有關(guān)。因此，每組數(shù)據(jù)集缺失率不同錯(cuò)誤率差別較大甚至成倍增加或減少。如表2所示，在Ronen數(shù)據(jù)集中，缺失率為1%和15%的均方殘值都為0.060，但是缺失率為1%的平均容積高缺失率15%近1倍，所以缺失率為1%的數(shù)據(jù)集錯(cuò)誤率小于缺失率為15%的數(shù)據(jù)集。

表2 四種實(shí)驗(yàn)數(shù)據(jù)集結(jié)果比較

實(shí)驗(yàn)表明， bi-SA方法采用k近鄰法、模擬退火法與雙聚類相結(jié)合，提高了數(shù)據(jù)填補(bǔ)的全局優(yōu)化能力。通過(guò)四組數(shù)據(jù)集實(shí)驗(yàn)，能夠取得較低的錯(cuò)誤率。在k近鄰法生成雙聚類種子的過(guò)程中，實(shí)驗(yàn)發(fā)現(xiàn)bi-SA方法對(duì)數(shù)據(jù)大小設(shè)定很敏感，對(duì)雙聚類容積和均方殘值影響很大，從而影響預(yù)估值的正確率。

4 結(jié) 語(yǔ)

本文提出的bi-SA方法是一種將k近鄰法、模擬退火法與雙聚類相結(jié)合填補(bǔ)算法，充分考慮了基因表達(dá)數(shù)據(jù)相關(guān)性特征，采用雙聚類質(zhì)量越高均方殘值越小填補(bǔ)的準(zhǔn)確性越高這一特性。通過(guò)與LLS、BPCA、Li-iLs和bi-BPCA四種方法的四組數(shù)據(jù)對(duì)比實(shí)驗(yàn)可知，bi-SA算法實(shí)驗(yàn)理想、填補(bǔ)準(zhǔn)確性較高。但是在實(shí)驗(yàn)過(guò)程中發(fā)現(xiàn)模擬退火法搜索過(guò)程的隨機(jī)性、初始溫度的設(shè)定和溫度降低變化的突跳性等缺陷，易導(dǎo)致實(shí)驗(yàn)結(jié)果產(chǎn)生大的波動(dòng)。因此，對(duì)提出方法還有待深入和完善，改進(jìn)的方法如對(duì)算法解的搜索性能進(jìn)一步實(shí)驗(yàn)、采用自適應(yīng)模擬退火法，根據(jù)雙聚類搜索進(jìn)展的反饋信息，自適應(yīng)地調(diào)整溫度變化和雙聚類的搜索強(qiáng)度等。

[1] Dudoit S,Yang Y H,Callow M J,et al.Statistical methods for identifying differentially expressed genes in replicated cDNA microarray experiments[J].Statistica Sinica,2002,12(1):111-140.

[2] 黃德雙.基因表達(dá)譜數(shù)據(jù)挖掘方法研究[M].北京:科學(xué)出版社,2009.

[3] Alizadeh A A,Eisen M B,Davis R E,et al.Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling[J].Nature,2000,403(6769):503-511.

[4] Bubin D B.Multiple imputation:a primer[J].Statistical Methods in Medical Research,1999,8(1):3-15.

[5] Rubin D B.Inference and missing data[J].Biometrika,1976,63(3):581-592.

[6] Jornsten R,Wang H Y,Welsh W J,et al.DNA microarray data imputation and significance analysis of differential expression[J].Bioinformatics,2005,21(22):4155-4161.

[7] Pan X,Tian Y,Huang Y,et al.Towards better accuracy for missing value estimation of epistatic miniarray profiling data by a novel ensemble approach[J].Genomics,2011,97(5):257-264.

[8] Wang Hsiuying,Chiu Chiachun,Wu Yiching,et al.Shrinkage regression-based methods for microarray missing value imputation[J].BMC Systems Biology,2013,7(6):16-18.

[9] Li H H,Shao F F,Li G Z.Semi-supervised imputation for microarray missing value estimation[C]//IEEE International Conference on Bioinformatics and Biomedicine.IEEE,2015:297-300.

[10] He C,Li H H,Zhao C,et al.Triple imputation for microarray missing value estimation[C]//IEEE International Conference on Bioinformatics and Biomedicine.IEEE,2015:208-213.

[11] Cheng Y,Church G M.Biclustering of Expression Data[C]//Proceedings of the Eighth International Conference on Intelligent Systems for Molecular Biology(ISMB),2000:93-103.

[12] Ji R,Liu Ding,Zhou Zuo.A Bicluster-Based Missing Value Imputation Method for Gene Expression Data[J].Journal of Computational Information Systems,2011,7(13):4810-4818.

[13] Meng Fangchi,Cheng Cai,Yan Hong.A bicluster-based bayesian principal component analysis method for microarray missing value estimation[J].IEEE journal of biomedical and health informatics,2014,18(3):863-871.

[14] Troyanskaya O,Cantor M,Sherlock G,et al.Missing value estimation methods for DNA microarrays[J].Bioinformatics,2001,17(6):520-525.

[15] 郝勝軒,宋宏,周曉鋒.一種基于雙聚類的缺失數(shù)據(jù)填補(bǔ)方法[J].計(jì)算機(jī)應(yīng)用研究,2015,32(3):674-678.

[16] Kirkpatrick S,Gelatt C D,Vecchi M P.Optimization by simulated annealing[J].Science,1983,220(4598):671-680.

[17] Baldwin D N,Vanchinathan V,Brown P O,et al.A gene expression program reflecting the innate immune response of cultured intestinal epithelial cells to infection by Listeria monocytogenes[J].Genome Biology,2003,4(1):1-14.

[18] Ronen M,Botstein D.Transcriptional Response of Steady-State Yeast Cultures to Transient Perturbations in Carbon Source[J].Proceedings of the National Academy of Sciences of the United States of America,2006,103(2):389-394.

[19] Ogawa N,Derisi J,Brown P O.New Components of a System for Phosphate Accumulation and Polyphosphate Metabolism in Saccharomyces cerevisiae Revealed by Genomic Expression Analysis[J].Molecular Biology of the Cell,2000,11(12):4309-4321.

[20] Yoshimoto H,Saltsman K,Gasch A P,et al.Genome-wide analysis of gene expression regulated by the calcineurin/Crz1p signaling pathway in Saccharomyces cerevisiae[J].Journal of Biological Chemistry,2002,277(34):31079-31088.

GENEDATAIMPUTATIONMETHODBASEDONSIMULATEANNEALINGOPTIMIZEDBICLUSTERING

Zhu Xian1Yang Ming1Ma Wei2Zhu Jun1

1(DepartmentofComputerScience,ZijinCollege,NanjingUniversityofScienceandTechnology,Nanjing210046,Jiangsu,China)2(SchoolofHotelManagement,NanjingInstituteofTourismandHospitality,Nanjing211100,Jiangsu,China)

Gene expression data are generated from the DNA microarray experiments of large data matrix, which can effectively extract the biological information. Due to the limitation of experimental conditions, gene expression data is often the presence of missing values and need to fill the missing data. The traditional missing data imputation method is based on the single feature of the gene expression data, and the correlation between the data matrix is not considered. The smaller the mean square residue of the dual clustering is, the higher the correlation between the data of the gene expression data is, and a method of missing data filling (bi-SA) based on simulated annealing optimized double clustering is proposed. This missing data imputation method used simulated annealing to optimize the optimal biclustering, and thus achieved the most effective imputation of missing data. Four groups of real gene expression data show that the bi-SA method has higher accuracy compared with other imputation methods.

Gene expression data Missing data imputation Simulated annealing algorithm Biclustering

2016-12-18。江蘇省高校自然科學(xué)研究項(xiàng)目(15KJB520017，16KJB520019)。朱嫻，講師，主研領(lǐng)域：智能優(yōu)化，生物信息。楊明，教授。馬衛(wèi)，講師。朱俊，講師。

TP391.9

A

10.3969/j.issn.1000-386x.2017.11.046