基于模擬退火優化雙聚類的基因數據填補方法

朱 嫻 楊 明 馬 衛 朱 俊

1(南京理工大學紫金學院計算機學院 江蘇 南京 210046) 2(南京旅游職業學院酒店管理學院 江蘇 南京 211100)

基于模擬退火優化雙聚類的基因數據填補方法

朱 嫻1楊 明1馬 衛2朱 俊1

1(南京理工大學紫金學院計算機學院 江蘇 南京 210046)2(南京旅游職業學院酒店管理學院 江蘇 南京 211100)

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

基因表達數據 缺失數據填補 模擬退火法 雙聚類

0 引 言

基因芯片技術[1]，又稱為DNA微陣列技術，是一種高通量的分子生物學手段，能夠對大量的基因進行快速測量，同時獲得成千上萬個基因在不同條件下的表達水平。基因芯片工作原理[2]是對兩個互補DNA分子雜交信號進行檢測分析，獲取待測生物樣本的基因表達豐度，在生物學和醫學研究等方面獲得了廣泛的應用，為探索生命科學提供了強有力的工具。從基因芯片實驗產生的基因表達數據通常是一個大形式的矩陣，不可避免地受到技術、實驗條件的限制而產生數據缺失問題。基因表達數據分析包括聚類、分類、特征值提取等方法，一般分析需要輸入完整的數據矩陣，對于具有缺失值的不完備基因表達數據矩陣不能直接處理。最初的解決方法是重復實驗找出缺失點或者刪除包含缺失值的基因，但是這兩種方法實用價值不高。隨后采用0值或行均值[3]填補缺失值的方法雖然簡單，但是處理數據的方式過于粗糙、誤差率高，對后續實驗結果產生較大影響，在底層的基因表達數據挖掘分析中被證明是缺少有用價值的。

數據缺失問題在各種應用研究領域都會出現，相關的缺失數據的填補方法研究歷史比較悠久，針對簡單填補法的不足，提出多重填補法MI(multiple imputation)[4-5]。通過對數據的填補、分析、合并，將數據進行綜合分析從而統計推斷出數據。近幾年，基因表達數據缺失值填補方法研究陸續提出，主要集中在混合全局的局部填補方法和雙聚類填補方法。

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

1 雙聚類基本原理

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

aij=aiJ+aIj-aIJ

(1)

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

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

(2)

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

(3)

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

(4)

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

(5)

2 基于模擬退火的雙聚類缺失數據方法

2.1 k近鄰法

k近鄰法是基因表達數據缺失值填補最常用方法之一，它是尋找與目標基因距離最近的k個具有完全值的基因。本文采用歐氏距離計算基因之間的距離，表示為：

(6)

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

2.2 模擬退火法

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

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

(7)

當ΔE<0，則接受種子S′作為當前新的雙聚類，否則計算狀態接受函數。

f=exp(-ΔE/T)

(8)

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

2.3 bi-SA方法描述

bi-SA步驟如下所述。

1) 找出基因表達數據集中缺失值x，計算x所在行的行均值，作為進行x初始值。

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

3) 模擬退火法優化雙聚類。

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

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

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

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

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

f) 逐漸減少溫度，然后轉b步。

4) 利用式(5)計算最優雙聚類缺失值x。

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

3 實驗結果與分析

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

表1 四種實驗數據集

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

(9)

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

圖1 不同數據集比較

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

表2 四種實驗數據集結果比較

實驗表明， bi-SA方法采用k近鄰法、模擬退火法與雙聚類相結合，提高了數據填補的全局優化能力。通過四組數據集實驗，能夠取得較低的錯誤率。在k近鄰法生成雙聚類種子的過程中，實驗發現bi-SA方法對數據大小設定很敏感，對雙聚類容積和均方殘值影響很大，從而影響預估值的正確率。

4 結 語

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

[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] 黃德雙.基因表達譜數據挖掘方法研究[M].北京:科學出版社,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] 郝勝軒,宋宏,周曉鋒.一種基于雙聚類的缺失數據填補方法[J].計算機應用研究,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。江蘇省高校自然科學研究項目(15KJB520017，16KJB520019)。朱嫻，講師，主研領域：智能優化，生物信息。楊明，教授。馬衛，講師。朱俊，講師。

TP391.9

A

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