基于改進標記傳播算法的基因表達譜數據分析

王年，葛芳，王俊生，唐俊

(安徽大學 計算智能與信號處理教育部重點實驗室，安徽 合肥，230039)

DNA 微陣列技術為腫瘤學提供了一種全新的研究手段，1 次基因芯片實驗可以獲得成千上萬個基因的信息。隨著DNA 微陣列技術的進步和儀器設備的更新，基因表達譜數據將不斷積累，這類數據的主要特點是樣本少、維數高、冗余基因和噪聲多，因此，如何對這類“新型”數據進行分析，挖掘其中包含的有效信息，已成為生物信息學研究的重點問題之一[1-3]。傳統的基因表達譜分析方法主要是對基因數據進行信息基因選取或特征屬性提取后，再用相應的分類或聚類方法對樣本進行識別[4-6]，其中聚類是一種無監督學習方法。傳統的聚類方法，如模糊C 均值[7]和自組織映射[8]等，均要求數據呈球狀分布，且這類算法易陷入局部最優。近年來，半監督聚類方法成為機器學習領域的研究熱點。Zhu 等[9]用高斯隨機場的方法處理半監督學習問題，提出了標記傳播算法，并根據隨機游走對標記傳播算法進行了概率解釋。Wang等[10-11]基于任意數據點的類標號可以由其近鄰點的類標號進行線性重構的假設，提出了基于線性近鄰的標記傳播算法；Bai 等[12]將標記傳播進行了擴展，提出了基于上下文分析的圖傳導算法，并將其應用于圖像分割、形狀檢索和匹配等領域。標記傳播是一種基于圖的半監督聚類方法，該方法利用少量已標記類別的樣本，通過傳播標記的方式來識別未知類別的樣本，能在任意形狀的樣本空間中進行聚類，克服了傳統聚類方法易陷入局部最優的缺點。然而，原始標記傳播算法迭代次數過大，每次迭代前都要重新標記已知類別的樣本，且最終的劃分準則也并不明確。針對上述問題，結合Gauss-Seidel 迭代[13]算法的思想，本文提出一種改進的標記傳播算法，并將其應用于基因表達譜數據分析。通過對白血病和結腸癌數據集的實驗，證明本文方法的有效性。

1 基因表達譜數據的權值矩陣表示

設{ g1, g2, …, gm}表示一個樣本中所有基因構成的集合，其中 gj(1≤ j ≤ m)表示1 個基因，m 表示基因個數；設{ x1, x2, …, xn}表示由所有實驗樣本構成的集合，其中 xi(1≤ i ≤ n)表示在同一條件下所有基因的表達值，n 表示樣本個數。由所有樣本構成的基因表達譜矩陣M 可表示為

其中：aij為基因gj在樣本xi中的表達值。

根據譜圖理論，以樣本為節點構造賦權圖G(V,E)(其中：V 為待聚類的樣本集；E=V×V 表示邊集)。設節點xi和xj之間的邊權值為Wij，本文所選擇的相似性度量為

其中： i ,j=1,2, … ,n；dij表示節點xi和xj之間的歐氏距離；σ 為權重參數，通常，σ 可以根據節點xi和xj的K 個近鄰的平均距離來自適應調整[14]，K 是一個經驗值。顯然，對于任意的xi，xj和σ ，有0≤Wij≤1，反映了節點xi和xj之間的親近程度，Wij越大，則節點xi和xj越有可能屬于同一類。另外，定義如下Laplace 矩陣：

其中：D 為度對角矩陣，且有 Dii=∑kWik。

2 標記傳播算法

標記傳播算法將某些樣本標記為已知，并將樣本集劃分為已標記樣本{ x1, x2, …, xl}和未標記樣本{xl+1,xl+2, …,xl+u}，其中，l + u=n，且滿足l＜＜ u。令YL=(y1, y2, …, yl)表示已知樣本的標記，同時定義一個標記序列f=[fLfU]T表示最終的標記結果，其中， fL=YL。標記傳播的步驟如下。

Step 1: 傳播標記，f←D-1Wf ；

Step 2: 強化標記數據， fL=YL；

Step 3: 返回Step 1，直到f 收斂。

通過step 1，將已標記樣本節點的標記傳播至相鄰的節點；同時，在標記傳播的過程中，已知樣本的標記是保持不變的，以保證標記數據的強度，指導標記傳播過程的正確性；算法迭代終止后，標記序列f包含了樣本的類別信息。Zhu[9]給出了標記傳播算法的收斂解：

其中：PUU和PUL是迭代矩陣P=D-1W 的子矩陣塊，即

因此，標記傳播算法的解將收斂于1 個固定解。需要指出的是，根據式(5)可以直接求出未知樣本的標記。

3 改進的標記傳播算法

通過上述方法任意數據點xi都可以映射為相應的實數值fi，但該算法在每次迭代前都要重新標記已知樣本，且由于標記傳播過程和標記強化過程分離，造成算法迭代次數過大；同時，該算法采用0-1 標記，需要選取適合的閾值對樣本進行分類。Zhu 等[9]使用的0.5 閾值缺乏穩定性，應用于復雜數據時并不能取得很好的效果。

3.1 基于Gauss-Seidel 迭代的標記傳播算法

設線性方程組Ax=b，矩陣A 為實對稱矩陣，且滿足A=D - W1-W2。其中：D 為對角矩陣；W1為下三角矩陣；W2為上三角矩陣。若迭代矩陣J=(D -W1)-1W2的譜半徑 ρ(J)＜1，則可以利用Gauss-Seidel迭代求解該方程組，如在第t+1次迭代后，該方程組的近似解為

令 x=f ,b=f0，其中： f0=(y1, …, yl,0, … ,0)，同時令A=IL+μL，μ∈(0,1)為可調參數，IL為對角矩陣，且滿足

根據Gauss-Seidel 迭代算法，經過t+1 次迭代后，數據點xi的標記為

由 于Aij=ILij+μ(Dij-Wij)=-μWij， 同 時 有Aii=ILii+μ(Dii-Wii)=ILii+μ Dii，因此，

由式(9)可知：在每次迭代過程中，數據點從初始標記節點中獲取一部分標記信息；同時，由于相似節點間的權重較大，使得數據點可以從其近鄰點中獲取一部分標記信息，當傳播終止后，相似節點的分布情況也趨于一致。化簡式(9)可得

式(10)可以改寫為

利用式(11)更新每個數據點的標記，直到收斂。這里，對初始標記點使用正負標記的方式，如對于兩類問題，將其中一類的若干個樣本標記為1，而另一類的若干個樣本標記為-1。對于最終的收斂結果，能以零為分割點對未知類別的樣本進行劃分。

3.2 收斂性證明

本節將證明標記序列f 的收斂性，展開式(11)可得

令Q=(IL+μD -μW1)-1μW2，顯然Qij≥ 0且∑jQij＜1。根據Perron-Frobenius 定理[13]，Q 的譜半徑小于1；同時，由于0＜μ＜1，且0＜Wij＜1，因此，

其中：E 為單位矩陣。因此， ft將最終收斂于

化簡式(15)，可得

因此，本文算法的解是收斂的。同樣，可以通過式(16)直接求出所有樣本的最終標記值。

4 實驗結果與分析

4.1 實驗數據

選用 2 組公開的癌癥數據集：白血病數據(leukemia)，共52 個樣本，其中24 個樣本為急性淋巴白血病(ALL)，28 個樣本為急性粒細胞白血病(AML)，每個樣本含有12 564 個基因表達數據(數據源于http://www.broad.mit.edu/cgi-bin/cancer/datasets.cgi) ；結腸癌數據(colon cancer)，共62 個樣本，其中22 個樣本被確定為正常樣本，40 個被確定為腫瘤樣本，每個樣本含有 2 000 個基因表達數據(數據源于http://linus.nci.nih.gov/～brb/DataArchive_New.html)。本文實驗是在酷睿雙核主頻2.60 GHz，內存為2 G 的計算機上運行的。

癌癥基因表達譜數據的獲取過程十分復雜，所得到的數據含有大量的噪聲，同時每個樣本都記錄了組織細胞中所有可測基因的表達水平，但只有較少數基因包含與類別相關的信息，因此，對數據進行預處理是必要的，定義下式對基因表達譜數據進行篩選：

其中：max(i)表示為第i 個基因在所有樣本中表達水平的最大值；min(i)表示為第i 個基因在所有樣本中表達水平的最小值；mean(i)表示為第i 個基因在所有樣本中表達水平的平均值；T 為給定的閾值。若某個基因的表達情況符合式(17)，則將該基因剔除。

4.2 白血病數據實驗

將白血病數據的第1 號樣本(ALL)標記為1，第25 和52 號樣本(AML)標記為-1，分別進行t=5，10，50，100 次迭代，更新每個樣本點的標記，觀察標記序列各分量的變化，結果如圖1 所示。

圖1 中，標記序列的前24 個分量對應的是白血病數據集的ALL 樣本，后28 個分量對應AML 樣本。由圖1 可知：當進行5 次迭代后，ALL 樣本的標記均大于0，而AML 樣本的標記均小于0，所有樣本被正確地劃分為2 類，證實本文方法可以快速且有效的實現樣本類別的劃分；隨著迭代次數的增加，2 類樣本標記值的區別更加明顯，同時第1 號樣本的標記值始終為1，而第25 和52 號樣本的標記值始終為-1，這說明初始標記點的標記信息在每次迭代時都被保留下來；當進行50 次和100 迭代后，標記序列各分量的變化并不明顯，說明經過較少次數的迭代就可以得到最終的收斂結果。

4.3 結腸癌數據實驗

將結腸癌數據的第1 號樣本(正常組織樣本)標記為1，第62 號樣本(腫瘤組織樣本)標記為-1，同樣分別進行t=5，20，100 和200 次迭代，結果如圖2 所示。

圖2 中，標記序列的前22 個分量對應的是結腸癌數據集的正常組織樣本，后40 個分量對應腫瘤組織樣本。由圖2 可知：初始標記點的標記值在標記傳播過程中始終保持不變。當經過5 次迭代后，2 類樣本標記的區別并不明顯，但經過20 次迭代后，就可以看出2 類樣本的標記的差異性，在正常組織樣本中，第18和20 號樣本被錯誤標記，而在腫瘤組織樣本中，第52，55 和58 號樣本被錯誤標記，準確率達到91.94%，在實際情況中，這些樣本中含有較多的偏離值和異常點，從而導致樣本異常表達；同時，進行100 次和200次迭代后，樣本標記值的變化已并不明顯，這同樣說明本文方法可以快速實現標記序列的收斂。

在理想情況下，初始標記點的選取是任意的，但在實際情況中，某些樣本含有較多的噪聲和異常表達值，若將其作為初始迭代點，則最終的迭代結果會受到一定的影響。因此，為了觀察初始標記點對最終迭代結果的影響，分別將結腸癌數據中被錯誤標記的第20，52 和55 號樣本標記為已知，結果如圖3 所示。

圖1 白血病數據的標記傳播過程Fig.1 Label propagation processes of leukemia data

圖2 結腸癌數據的標記傳播過程Fig.2 Label propagation processes of colon cancer data

圖3 異常樣本被標記為已知的迭代結果Fig.3 Iterative results by labeling abnormal samples

圖3(a)中僅標記了第1 和23 號樣本；圖3(b)中除標記第1 和23 號樣本外，同時將第20 號樣本標記為1；圖3(c)和(d)除標記第1 和23 號樣本外，分別將第52 號和55 號樣本標記為-1。由圖3(b)可知：當將第20 號樣本標記為已知，除20 號樣本外的正常組織樣本的標記值并沒有受到明顯的影響，腫瘤組織樣本中的第52，55 和58 號樣本的標記值也基本沒有變化，而其他樣本標記值的絕對值減小卻趨于減小，同時第24 和59 號樣本被錯誤標記到另一類中。由圖3(c)～(d)可知：當將第52 或55 號樣本標記為已知時，第13至22 號樣本全部被錯誤標記，同時第1 至12 號樣本的標記值的絕對值也相應地減小，這主要是由于在標記傳播過程中，某個樣本的標記是根據其近鄰樣本點的標記進行更新的。因此，若將某個異常表達的樣本標記為已知，則與該樣本相似的正常表達樣本，其標記值也會趨向于與該樣本的相似，甚至會被錯誤標記。

4.4 對比驗證

為了進一步驗證本文方法的正確性，將本文方法與傳統的S2N_KNN[1]法以及另外2 種基于圖的方法即局部線性嵌入法(Locally Linear Embedding,LLE)[5,15]和歸一化割法(Normalized Cut, Ncut)[16]進行對比。S2N_KNN 以“信噪比”為指標選取特征基因，再用K-近鄰(KNN)分類器對樣本進行分類；LLE 是一種流行學習方法，本文先以LLE 提取樣本特征，再結合KNN 實現樣本的分類。Ncut 通過計算規范Laplace矩陣的次小特征值對應的特征向量(Fiedler 向量)來進行聚類，結果如表1 所示。

表1 本文方法與另外3 種方法的實驗結果比較Table 1 Comparisons of experiment results with other three methods

由表1 可知：無論從準確率還是運算效率方面，本文方法都比傳統的S2N_KNN 法的高。傳統方法以降維來提取特征基因勢必會丟失部分含有分類信息的基因，而本文方法將樣本數據作為高維空間中的點，所構造的權值矩陣反映了樣本之間的關系，使得原來的樣本分類信息完全映射到低維的權值矩陣中。同時，本文方法首先以樣本為節點構圖，其低樣本特性決定了構造的矩陣規模較小，從而具有較低的運算復雜度，有利于基因表達數據的快速聚類。

另外，從準確率方面，本文方法明顯高于另外2種基于圖的方法。LLE 算法是一種有效的可視化方法，但其要求數據服從流形分布，提取的特征并不具有很強的分類能力。Ncut 在聚類時傾向于所得到的類具有相同的類內相似度，且Fiedler 向量并不一定能夠反映原始數據的結構；而本文方法通過在近鄰點之間傳播標記的方式來學習分類，并不會受到數據分布形狀的限制。同時，在效率上，本文方法并不具有明顯的優勢，這是因為這幾種方法都是利用構圖的方法將高維樣本映射到低維空間，構造的矩陣規模是相同的。

5 結論

(1) 針對原始標記傳播算法存在的問題，結合Gauss-Seidel 迭代的思想，提出了一種改進的標記傳播算法，并將其應用于基因表達譜數據分析。

(2) 通過2 組癌癥數據集的實驗，證明本文方法可以快速有效地實現基因表達譜數據的聚類。該方法克服了迭代次數過大和重復標記數據點的缺點；同時在數據類別劃分時使用正負標記的方式，避免了采用0-1 標記時閾值選取的不確定性。

(3) 與傳統標記傳播算法一樣，該方法的不足之處是初始標記點的選取對最終迭代結果具有一定的影響，這有待進一步研究與探討。

[1] Singh D, Febbo P G, Ross K, et al. Gene expression correlates of clinical prostate cancer behavior[J]. Cancer Cell, 2002, 1(2):203-209.

[2] Dash S, Patra B. A study on gene selection and classification algorithms for classification of gene expression profile[J].International Journal of Research and Reviews in Computer Science, 2011, 2(5): 1212-1217.

[3] 沈威, 鄭明, 劉桂霞, 等. 基于奇異值求通解方法進行基因調控網絡構建[J]. 中南大學學報(自然科學版), 2012, 43(4):1377-1381.SHEN Wei, ZHENG Ming, LIU Guixia, et al. Construction of gene regulatory network based on singular value decomposition for solution set[J]. Journal of Central South University (Science and Technology), 2012, 43(4): 1377-1381.

[4] 李宏, 李翔, 吳敏, 等. 基于閉合模式的高維基因表達譜多類分類[J]. 中南大學學報(自然科學版), 2008, 39(5): 1035-1041.LI Hong, LI Xiang, WU Min, et al. Multi-class classification of high-dimension gene expression profile based on closed patterns[J]. Journal of Central South University (Science and Technology), 2008, 39(5): 1035-1041.

[5] LI Bo, ZHENG Chunhou, HYANG Deshuang, et al. Gene expression data classification using locally linear discriminant embedding[J]. Computers in Biology and Medicine, 2010,40(10): 802-810.

[6] Kancherla K, Mukkamala S. Feature selection for lung cancer detection using SVM based recursive feature elimination method[J]. Machine Learning and Data Mining in Bioinformatics, 2012, 7246: 168-176.

[7] Tari L, Baral C, Kim S. Fuzzy c-means clustering with prior biological knowledge[J]. Journal of Biomedical Informatics,2009, 42(1): 74-81.

[8] Patterson A D, LI Henghong, Eichler G S, et al.UPLC-ESI-TOFMS-based metabolomics and gene expression dynamics inspector self-organizing metabolomic maps as tools for understanding the cellular response to ionizing radiation[J].American Chemical Society, 2008, 80(3): 665-674.

[9] ZHU Xiaojin. Semi-supervised learning with graphs[D].Pennsylvania: Carnegie Mellon University. School of Computer Science, 2005: 5-8.

[10] WANG Fei, ZHANG Changshui. Label propagation through linear neighborhoods[J]. IEEE Transactions on Knowledge and Data Engineering, 2008, 20(1): 55-67.

[11] WANG Jingdong, WANG Fei, ZHANG Changshui, et al. Linear neighborhood propagation and its applications[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2009,31(9): 1600-1615.

[12] BAI Xiang, YANG Xingwei, Latecki L J, et al. Learning context-sensitive shape similarity by graph transduction[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2010,32(5): 861-874.

[13] Saad Y. Iterative methods for sparse linear systems[M]. Boston:PWS Publishing Company, 1996: 104-107.

[14] Zelnik-Manor L, Perona P. Self-tuning spectral clustering[J].Advances in Neural Information Processing Systems, 2004,17(2): 1601-1608.

[15] Shi C, Chen L. Feature dimension reduction for microarray data analysis using locally linear embedding[C]// Proceeding of 3rd Asia-Pacific Bioinformatics Conference. London: Imperial College Press, 2005: 211-217.

[16] Dhillon I S, GUAN Yuqiang, Kulis B. Kernel k-means-spectral clustering and normalized cuts[C]// Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. New York: Association for Computing Machinery, 2004: 551-556.