采用自由能方法預測RNA二級結構時, 如何精確有效地從次優結構中篩選出真實的二級結構成為RNA結構預測中的關鍵。采用聚類技術對次優結構集合進行分析, 可有效地提高預測結果的精度。本文利用RBP分數矩陣, 提出一種基于增量中心候選集的改進k-medoids算法。它將隨機選擇初始中心并進行首次劃分后以中心候選集逐一擴展的方式進行中心輪換, 以降低算法的復雜度。實驗表明, 該算法能取得更高的CH值, 且能有效地縮短計算時間。
引用本文： 王常武, 劉小鳳, 王寶文, 劉文遠. IC-kmedoids:適用于RNA二級結構預測的聚類算法. 生物醫學工程學雜志, 2015, 32(1): 99-103. doi: 10.7507/1001-5515.20150018 復制
引言
RNA是組成生命體最重要的高分子化合物之一，它可以通過一些特定功能影響基因的表達，控制蛋白質合成，從而影響生命體的生長發育[1]。尤其是人類免疫缺陷病毒(human immunodeficiency virus，HIV)、肝炎病毒等，其病毒功能直接由RNA決定，而非DNA。RNA分子的功能特性與其結構緊密相關[2]，所以科研人員更著重于研究其結構特性。現階段對于RNA結構研究的關鍵點是對其二級結構的研究。
RNA二級結構預測算法大體上可以分為兩類：第一類是基于比較序列分析的方法，第二類是基于最小自由能的方法。基于序列比對的方法原理是先對多條序列中的互補堿基進行配對，然后在已知的數據庫中查找與被測序列結構相似的序列，從而來推斷未知序列的二級結構[3]。該方法適用于具有一系列同源RNA分子的情況。然而，當沒有任何先驗知識時，預測RNA二級結構一般采用最小自由能模型[4]。該模型假設真實的RNA會折疊成一個具有最小自由能的結構，然后利用動態規劃算法獲取此結構。然而實驗證明，RNA真實結構往往不是具有最小自由能的結構，為解決真實結構與最小自由能結構不一致的問題，Zuker等[5]提出次優結構的概念。他認為，真實二級結構的自由能也許不是最小的，但也應該具有一個較小的值使得RNA分子相對穩定。因此，可以人為設定一個閾值，與自由能相差該閾值以內的所有二級結構都有可能是真實結構。然后，將次優結構全部交由生物人員進行鑒定。Ding等[6]提出對次優結構集合進行聚類處理的方法，然后選出每一簇的代表結構，進行下一步實驗，以提高工作效率。他首次提出采用BP距離，即不一致配對的數目作為兩結構間的距離。計算兩兩結構間的堿基對距離構成堿基對矩陣(此矩陣為對稱矩陣)，以此矩陣作為輸入，進行裂解分層聚類。Agius等[7]在此基礎上又提出改進的堿基對距離，即松馳堿基對(relaxed base pair, RBP)距離。求出RBP矩陣，對其進行k-均值聚類，取得了比文獻[7]較優的效果。
最小自由能方法中聚類技術的引入，使得RNA二級結構預測的精度和效率都有了很大的提高，為研究者提供了一種新的思路。然而，針對RNA預測的聚類算法還很少。層次聚類中的數據通過算法構造樹連接，適用于比較同源或是相似的RNA結構。k-均值算法對異常數據很敏感，而且初始中心的選取對聚類結果有很大的影響。所以，有效的聚類算法對于RNA結構預測至關重要。由于現階段科研水平的限制，研究人員對RNA二級結構的認識還不是很清楚，且生物數據量龐大，所以在實驗的過程中很難避免異常數據。傳統的k-medoids算法能有效地處理異常數據，但計算量偏大。
針對這個問題，本文提出了一種基于增量中心候選集的改進k-medoids算法—IC-kmedoids。它將隨機選擇初始中心并進行首次劃分后，以中心候選集逐一擴大的方式進行中心輪換，以降低算法的復雜度。該算法不僅延續了原始k-medoids算法的優勢，保證了聚類效果的優越性與穩定性，同時算法效率也取得了較大程度的提高。采用CH(Calinski and Harabasz)索引[8]在[2,14]區間中尋找最佳聚類數目k。對于選取N個結構的序列，構造N×N矩陣，計算兩兩結構間的RBP距離，構成矩陣中的元素(矩陣為對稱矩陣)。矩陣中的每行作為一個數據，每一列作為數據的特征，對這N個數據進行聚類分析。實驗證明，IC-kmedoids算法能快速、有效地將RNA折疊結構聚類。
1 基于增量中心候選集的改進k-medoids算法
1.1 采用基于增量中心候選集的中心替換策略
對于某一中心點，若要對其進行調整，一些較優的候選中心點最有可能是該中心附近區域中的點，如本簇或其附近簇中的點，因此可以在搜索新的中心點時限定搜索范圍[9]。具體方法為：假設k為目標聚類簇數，則中心替換分為k次進行，首先在每次替換前生成中心候選集，且中心候選集在迭代過程中不斷增大。假設當前替換為第i次，則對所有中心點進行替換時，設定其搜索范圍為離該中心點最近的i個簇(包含被替換中心所在簇)中的所有非中心對象，從而保證新中心點候選集從1, 2，…，k個子集不斷增大，直至最后的搜索范圍擴大到所有的非中心點，與原始k-medoids算法保持一致。其算法描述如下：
輸入：包含N個對象的數據集合D，以及結果簇的數目k。
輸出：k個簇，使得所有對象與所屬簇的中心點距離總和最小。
(1)在N個數據中隨意選擇k個數據作為初始的中心點集{O1, O2, …, Ok}；
(2)以當前中心點集合{O1, O2, …, Ok}對其余對象進行初始劃分；
(3)執行增量中心候選集替換策略
令i=1；
while i≤k
(輪換中心并計算代價；代價計算：新中心與本簇中所有數據點距離平方和-原始數據O與本簇所有數據點距離平方和)
repeat
??選擇一個未選擇過的中心點O；
??選出距O最近的i個簇(包含O本身所在簇)內的所有非中心對象，生成替換O的候選中心集合C；
??repeat
??在C中選擇一個未選擇過的數據點d，計算用d代替O的代價并記錄在S中；
??until C中的數據對象都被選擇過；
?until中心點集合{O1, O2, …, OK}中的所以數據都被選擇過；
(實施中心替換)
?if min(S)＜0(min(S)表示代價集合S中的最小值)
?用min(S)對應的數據對象與中心點進行替換，形成一個新的k個中心點的集合{O1, O2, …, Ok}；
?把剩余的n-k個對象指派給離它最近的中心點所代表的簇；
?返回while處重新開始執行；
else(即min(S)≥0時, 不進行替換)
?i=i+1
end if
end while
(4) 算法結束
1.2 聚類數目k值的產生——CH索引
本文提出的算法要產生特定數目的簇，然而算法本身并不產生聚類數目k。CH提供一種度量方式：一個好的聚類算法應使得類內距離相近，且類間距離增大。本文中樣本間的距離采用歐氏距離進行計算。假定N個實例被分成k類，，其中B(k)為類間距離平方和，W(k)為類內距離平方和。B(k)=(EC, CCi)，W(k)=(CCi, Iij)，其中ni為第i簇中二級結構的數目，D(EC, CCi)為整體的聚類中心(EC)與第i簇的聚類中心(CCi)之間的距離，D(CCi, Iij)為簇i中的第j個數據點與本簇中心點之間的距離。k值默認為1，由于CH索引并沒有全局最大的值，因此在k≥2的區間內，進行搜索驗證，以獲得最大CH。本文中k的取值范圍定為[2,14]。
2 RBP矩陣說明及其算法
假設一條RNA序列長度為n，第i個堿基和第j個堿基分別表示為ri和rj，其中i＜j。采用螺旋區點陣圖的方式來描繪RNA二級結構。構造點陣圖，i(5端)從上到下依次增長，j(3端)從左到右依次增長。給定該序列的兩個結構S1和S2，S1中的堿基配對采用紅點表示，S2中的堿基配對采用綠點表示。當一個配對同時出現在S1和S2中時，用黑點表示。堿基對標準是統計紅綠點總數，使用一套嚴格的相近性標準。RBP引入一個松弛變量t，S1和S2兩結構間的RBP距離記為ρt(S1, S2)。若兩個結構是相同的，則所有點都是黑色的，RBP=0。否則，挑選一個與其不同色的距離最遠的紅色或是綠色點進行排除。之后，在每個配對位置放置一個t×t方格(紅色點位置放紅色方格，綠色點位置放綠色方格，黑色點位置放黑色方格)。不考慮被排除的配對點，如果每一個紅綠小格都與其它顏色小格有重疊(包括邊或是頂點)，則RBP距離設為1；如果，還有未重疊的，那么排除它，并把方格擴大到2t×2t；循環上述過程至d，直到紅綠小格都與其它顏色小格有重疊時結束，小方格大小為td×td，此時，RBP為d。
若t=0，這些小方格呈現為單獨的一點，因此所有的紅色和綠色的點都會被排除，最終只剩下黑色共有的堿基配對。這表明，在t=0時，RBP標準與堿基對標準是相同的。隨著t的增加，RBP距離保持不變或是減小。如果，t的大小設置為序列的長度，若兩結構相同，則沒有配對堿基被排除，RBP=0；否則，在排除一個配對方格后，小方格變為t×t，紅綠小方格都相互交叉，RBP=1。當t足夠大時，RBP變得不再重要，相同結構的RBP距離為0，不同結構的為1。
為了說明問題，本文選取一個較簡單序列的兩個折疊結構為例，如圖 1所示。在兩個結構中同時出現的配對用黑色標記，只出現在左圖結構中的配對采用紅色標記，只出現在右圖結構中的配對用綠色標記，未配對堿基用灰色標記。以t=1為例，圖 2為該序列兩結構的點陣圖，方格大小為2×2，所以RBP距離為2。
圖1
同一序列的兩個折疊結構
Figure1.
Two secondary structures for the same RNA sequence
RBP算法描述如下：假設i·j表示核苷酸i和j配對，則設兩對配對堿基i·j與i′·j′之間的距離為δbb(i·j, i′·j′)=max{｜i-i′｜, ｜j-j′｜}。例G23C94和U30A86的配對距離為8，即max{30-23, 94-86}。一對配對堿基i·j與另一結構S之間的距離δbs定義為δbs(i·j, S)=mini′·j′∈S{δbb(i·j, i′·j′)}。假設兩個二級結構S1和S2，分別包含M1和M2對配對堿基。則兩結構間所有配對距離為：S1中M1個堿基配對距S2結構的距離，與S2中M2個堿基配對距S1結構的距離。這些距離并不是完全不同的。若S1結構和S2結構中含有M3對一致的堿基配對時，那么會出現M3個0。對這些距離進行降序排列操作，記作{Δ1, Δ2, …, ΔM1+M2}，即對于任意i＜j，則有Δi≥Δj。對于任意t≥0，定義ρt(S1, S2)=min{m∈Z｜m≥0, Δk≤tm若k＞m}。當t=0時，tm恒為0，此時RBP標準等同于堿基對標準，即ρ0(S1, S2)為距離不為0的堿基配對數目。如果t＞maxΔi, 則m為1，此時，。可以設置不同的t值來計算RBP距離。本文中選取的值為1。
3 實現方法
3.1 數據集的選擇
從NCBI(美國國立生物技術中心)的Reference Sequence(RefSeq)數據庫[10]中挑選4條無冗余的人源mRNA序列(NM_130776、NM_016815、NM_013961、NM_002231)，兩條果蠅mRNA序列(NM_143005、NM_206557)和一條褐家鼠mRNA序列(NM_031515)，見表 1。對于每一條序列，從次優結構集合中隨機抽取1 000個折疊結構進行聚類。
表1
實驗數據
Table1.
Experimental data set
3.2 實驗步驟
(1)對于每一條RNA序列的1 000個折疊結構集合{S1, S2, …, S1000}，利用RBP算法求出第i個(1≤i≤1 000)結構與所有結構之間的RBP距離，作為此結構的屬性特征，即第i行第j列(1≤j≤1 000)的元素為ρ1(Si, Sj)，利用此1 000×1 000矩陣作為數據集D(主對角線上元素為零，相同結構的RBP距離為零)。
(2)在k∈[2,14]的區間上，分別采用改進的k-medoids算法對數據集D進行聚類，在聚類的過程中利用歐幾里得距離衡量數據之間的相似性。并求取各自所對應的CH值，使得CH值取得最大值時的k設定為聚類的最終數目K。圖 3為NM_130776 XAGE3的CH索引圖，其最終結果簇數設定為3。
圖3
NM_130776 XAGE3的CH索引圖
Figure3.
CH index of NM_130776 XAGE3
(3)采用主成分分析方法分析數據集D，挑選出第一主成分和第二主成分構成數據集合D′。利用改進的k-medoids算法把數據集D′分為K類，并用圖形可視化。圖 4為NM_130776 XAGE3的聚類結果圖。
圖4
NM_130776 XAGE3聚類結果圖
Figure4.
Results of clustering on NM_130776 XAGE3
3.3 實驗結果比對分析
采用原k-medoids算法與本文提出的算法進行比較，實驗結果見表 2。
表2
實驗結果
Table2.
Test results
由表 2可以看出，針對不同物種的mRNA序列，本文提出的算法取得的最大CH值均高于原始k-medoids算法的結果。這說明針對不同的物種，IC-kmedoids算法能在類間距離增大的同時，有效地減小類內距離。改進算法的運行時間比原始算法速率提高約30%。實驗結果表明，基于增量中心候選集的算法可快速、有效地對RNA折疊結構進行聚類。
4 結論
本文提出的基于增量中心候選集IC-kmedoids算法，能有效處理RNA折疊結構的聚類問題。該算法在中心點替換的過程中采取中心候選集逐一擴大的方式進行計算，有效地提高了算法的運行速率，并保證了聚類的效果。實驗表明，IC-kmedoids算法對于解決高維、復雜的RNA生物數據的聚類問題十分有效。
引言
RNA是組成生命體最重要的高分子化合物之一，它可以通過一些特定功能影響基因的表達，控制蛋白質合成，從而影響生命體的生長發育[1]。尤其是人類免疫缺陷病毒(human immunodeficiency virus，HIV)、肝炎病毒等，其病毒功能直接由RNA決定，而非DNA。RNA分子的功能特性與其結構緊密相關[2]，所以科研人員更著重于研究其結構特性。現階段對于RNA結構研究的關鍵點是對其二級結構的研究。
RNA二級結構預測算法大體上可以分為兩類：第一類是基于比較序列分析的方法，第二類是基于最小自由能的方法。基于序列比對的方法原理是先對多條序列中的互補堿基進行配對，然后在已知的數據庫中查找與被測序列結構相似的序列，從而來推斷未知序列的二級結構[3]。該方法適用于具有一系列同源RNA分子的情況。然而，當沒有任何先驗知識時，預測RNA二級結構一般采用最小自由能模型[4]。該模型假設真實的RNA會折疊成一個具有最小自由能的結構，然后利用動態規劃算法獲取此結構。然而實驗證明，RNA真實結構往往不是具有最小自由能的結構，為解決真實結構與最小自由能結構不一致的問題，Zuker等[5]提出次優結構的概念。他認為，真實二級結構的自由能也許不是最小的，但也應該具有一個較小的值使得RNA分子相對穩定。因此，可以人為設定一個閾值，與自由能相差該閾值以內的所有二級結構都有可能是真實結構。然后，將次優結構全部交由生物人員進行鑒定。Ding等[6]提出對次優結構集合進行聚類處理的方法，然后選出每一簇的代表結構，進行下一步實驗，以提高工作效率。他首次提出采用BP距離，即不一致配對的數目作為兩結構間的距離。計算兩兩結構間的堿基對距離構成堿基對矩陣(此矩陣為對稱矩陣)，以此矩陣作為輸入，進行裂解分層聚類。Agius等[7]在此基礎上又提出改進的堿基對距離，即松馳堿基對(relaxed base pair, RBP)距離。求出RBP矩陣，對其進行k-均值聚類，取得了比文獻[7]較優的效果。
最小自由能方法中聚類技術的引入，使得RNA二級結構預測的精度和效率都有了很大的提高，為研究者提供了一種新的思路。然而，針對RNA預測的聚類算法還很少。層次聚類中的數據通過算法構造樹連接，適用于比較同源或是相似的RNA結構。k-均值算法對異常數據很敏感，而且初始中心的選取對聚類結果有很大的影響。所以，有效的聚類算法對于RNA結構預測至關重要。由于現階段科研水平的限制，研究人員對RNA二級結構的認識還不是很清楚，且生物數據量龐大，所以在實驗的過程中很難避免異常數據。傳統的k-medoids算法能有效地處理異常數據，但計算量偏大。
針對這個問題，本文提出了一種基于增量中心候選集的改進k-medoids算法—IC-kmedoids。它將隨機選擇初始中心并進行首次劃分后，以中心候選集逐一擴大的方式進行中心輪換，以降低算法的復雜度。該算法不僅延續了原始k-medoids算法的優勢，保證了聚類效果的優越性與穩定性，同時算法效率也取得了較大程度的提高。采用CH(Calinski and Harabasz)索引[8]在[2,14]區間中尋找最佳聚類數目k。對于選取N個結構的序列，構造N×N矩陣，計算兩兩結構間的RBP距離，構成矩陣中的元素(矩陣為對稱矩陣)。矩陣中的每行作為一個數據，每一列作為數據的特征，對這N個數據進行聚類分析。實驗證明，IC-kmedoids算法能快速、有效地將RNA折疊結構聚類。
1 基于增量中心候選集的改進k-medoids算法
1.1 采用基于增量中心候選集的中心替換策略
對于某一中心點，若要對其進行調整，一些較優的候選中心點最有可能是該中心附近區域中的點，如本簇或其附近簇中的點，因此可以在搜索新的中心點時限定搜索范圍[9]。具體方法為：假設k為目標聚類簇數，則中心替換分為k次進行，首先在每次替換前生成中心候選集，且中心候選集在迭代過程中不斷增大。假設當前替換為第i次，則對所有中心點進行替換時，設定其搜索范圍為離該中心點最近的i個簇(包含被替換中心所在簇)中的所有非中心對象，從而保證新中心點候選集從1, 2，…，k個子集不斷增大，直至最后的搜索范圍擴大到所有的非中心點，與原始k-medoids算法保持一致。其算法描述如下：
輸入：包含N個對象的數據集合D，以及結果簇的數目k。
輸出：k個簇，使得所有對象與所屬簇的中心點距離總和最小。
(1)在N個數據中隨意選擇k個數據作為初始的中心點集{O1, O2, …, Ok}；
(2)以當前中心點集合{O1, O2, …, Ok}對其余對象進行初始劃分；
(3)執行增量中心候選集替換策略
令i=1；
while i≤k
(輪換中心并計算代價；代價計算：新中心與本簇中所有數據點距離平方和-原始數據O與本簇所有數據點距離平方和)
repeat
??選擇一個未選擇過的中心點O；
??選出距O最近的i個簇(包含O本身所在簇)內的所有非中心對象，生成替換O的候選中心集合C；
??repeat
??在C中選擇一個未選擇過的數據點d，計算用d代替O的代價并記錄在S中；
??until C中的數據對象都被選擇過；
?until中心點集合{O1, O2, …, OK}中的所以數據都被選擇過；
(實施中心替換)
?if min(S)＜0(min(S)表示代價集合S中的最小值)
?用min(S)對應的數據對象與中心點進行替換，形成一個新的k個中心點的集合{O1, O2, …, Ok}；
?把剩余的n-k個對象指派給離它最近的中心點所代表的簇；
?返回while處重新開始執行；
else(即min(S)≥0時, 不進行替換)
?i=i+1
end if
end while
(4) 算法結束
1.2 聚類數目k值的產生——CH索引
本文提出的算法要產生特定數目的簇，然而算法本身并不產生聚類數目k。CH提供一種度量方式：一個好的聚類算法應使得類內距離相近，且類間距離增大。本文中樣本間的距離采用歐氏距離進行計算。假定N個實例被分成k類，，其中B(k)為類間距離平方和，W(k)為類內距離平方和。B(k)=(EC, CCi)，W(k)=(CCi, Iij)，其中ni為第i簇中二級結構的數目，D(EC, CCi)為整體的聚類中心(EC)與第i簇的聚類中心(CCi)之間的距離，D(CCi, Iij)為簇i中的第j個數據點與本簇中心點之間的距離。k值默認為1，由于CH索引并沒有全局最大的值，因此在k≥2的區間內，進行搜索驗證，以獲得最大CH。本文中k的取值范圍定為[2,14]。
2 RBP矩陣說明及其算法
假設一條RNA序列長度為n，第i個堿基和第j個堿基分別表示為ri和rj，其中i＜j。采用螺旋區點陣圖的方式來描繪RNA二級結構。構造點陣圖，i(5端)從上到下依次增長，j(3端)從左到右依次增長。給定該序列的兩個結構S1和S2，S1中的堿基配對采用紅點表示，S2中的堿基配對采用綠點表示。當一個配對同時出現在S1和S2中時，用黑點表示。堿基對標準是統計紅綠點總數，使用一套嚴格的相近性標準。RBP引入一個松弛變量t，S1和S2兩結構間的RBP距離記為ρt(S1, S2)。若兩個結構是相同的，則所有點都是黑色的，RBP=0。否則，挑選一個與其不同色的距離最遠的紅色或是綠色點進行排除。之后，在每個配對位置放置一個t×t方格(紅色點位置放紅色方格，綠色點位置放綠色方格，黑色點位置放黑色方格)。不考慮被排除的配對點，如果每一個紅綠小格都與其它顏色小格有重疊(包括邊或是頂點)，則RBP距離設為1；如果，還有未重疊的，那么排除它，并把方格擴大到2t×2t；循環上述過程至d，直到紅綠小格都與其它顏色小格有重疊時結束，小方格大小為td×td，此時，RBP為d。
若t=0，這些小方格呈現為單獨的一點，因此所有的紅色和綠色的點都會被排除，最終只剩下黑色共有的堿基配對。這表明，在t=0時，RBP標準與堿基對標準是相同的。隨著t的增加，RBP距離保持不變或是減小。如果，t的大小設置為序列的長度，若兩結構相同，則沒有配對堿基被排除，RBP=0；否則，在排除一個配對方格后，小方格變為t×t，紅綠小方格都相互交叉，RBP=1。當t足夠大時，RBP變得不再重要，相同結構的RBP距離為0，不同結構的為1。
為了說明問題，本文選取一個較簡單序列的兩個折疊結構為例，如圖 1所示。在兩個結構中同時出現的配對用黑色標記，只出現在左圖結構中的配對采用紅色標記，只出現在右圖結構中的配對用綠色標記，未配對堿基用灰色標記。以t=1為例，圖 2為該序列兩結構的點陣圖，方格大小為2×2，所以RBP距離為2。
圖1
同一序列的兩個折疊結構
Figure1.
Two secondary structures for the same RNA sequence
RBP算法描述如下：假設i·j表示核苷酸i和j配對，則設兩對配對堿基i·j與i′·j′之間的距離為δbb(i·j, i′·j′)=max{｜i-i′｜, ｜j-j′｜}。例G23C94和U30A86的配對距離為8，即max{30-23, 94-86}。一對配對堿基i·j與另一結構S之間的距離δbs定義為δbs(i·j, S)=mini′·j′∈S{δbb(i·j, i′·j′)}。假設兩個二級結構S1和S2，分別包含M1和M2對配對堿基。則兩結構間所有配對距離為：S1中M1個堿基配對距S2結構的距離，與S2中M2個堿基配對距S1結構的距離。這些距離并不是完全不同的。若S1結構和S2結構中含有M3對一致的堿基配對時，那么會出現M3個0。對這些距離進行降序排列操作，記作{Δ1, Δ2, …, ΔM1+M2}，即對于任意i＜j，則有Δi≥Δj。對于任意t≥0，定義ρt(S1, S2)=min{m∈Z｜m≥0, Δk≤tm若k＞m}。當t=0時，tm恒為0，此時RBP標準等同于堿基對標準，即ρ0(S1, S2)為距離不為0的堿基配對數目。如果t＞maxΔi, 則m為1，此時，。可以設置不同的t值來計算RBP距離。本文中選取的值為1。
3 實現方法
3.1 數據集的選擇
從NCBI(美國國立生物技術中心)的Reference Sequence(RefSeq)數據庫[10]中挑選4條無冗余的人源mRNA序列(NM_130776、NM_016815、NM_013961、NM_002231)，兩條果蠅mRNA序列(NM_143005、NM_206557)和一條褐家鼠mRNA序列(NM_031515)，見表 1。對于每一條序列，從次優結構集合中隨機抽取1 000個折疊結構進行聚類。
表1
實驗數據
Table1.
Experimental data set
3.2 實驗步驟
(1)對于每一條RNA序列的1 000個折疊結構集合{S1, S2, …, S1000}，利用RBP算法求出第i個(1≤i≤1 000)結構與所有結構之間的RBP距離，作為此結構的屬性特征，即第i行第j列(1≤j≤1 000)的元素為ρ1(Si, Sj)，利用此1 000×1 000矩陣作為數據集D(主對角線上元素為零，相同結構的RBP距離為零)。
(2)在k∈[2,14]的區間上，分別采用改進的k-medoids算法對數據集D進行聚類，在聚類的過程中利用歐幾里得距離衡量數據之間的相似性。并求取各自所對應的CH值，使得CH值取得最大值時的k設定為聚類的最終數目K。圖 3為NM_130776 XAGE3的CH索引圖，其最終結果簇數設定為3。
圖3
NM_130776 XAGE3的CH索引圖
Figure3.
CH index of NM_130776 XAGE3
(3)采用主成分分析方法分析數據集D，挑選出第一主成分和第二主成分構成數據集合D′。利用改進的k-medoids算法把數據集D′分為K類，并用圖形可視化。圖 4為NM_130776 XAGE3的聚類結果圖。
圖4
NM_130776 XAGE3聚類結果圖
Figure4.
Results of clustering on NM_130776 XAGE3
3.3 實驗結果比對分析
采用原k-medoids算法與本文提出的算法進行比較，實驗結果見表 2。
表2
實驗結果
Table2.
Test results
由表 2可以看出，針對不同物種的mRNA序列，本文提出的算法取得的最大CH值均高于原始k-medoids算法的結果。這說明針對不同的物種，IC-kmedoids算法能在類間距離增大的同時，有效地減小類內距離。改進算法的運行時間比原始算法速率提高約30%。實驗結果表明，基于增量中心候選集的算法可快速、有效地對RNA折疊結構進行聚類。
4 結論
本文提出的基于增量中心候選集IC-kmedoids算法，能有效處理RNA折疊結構的聚類問題。該算法在中心點替換的過程中采取中心候選集逐一擴大的方式進行計算，有效地提高了算法的運行速率，并保證了聚類的效果。實驗表明，IC-kmedoids算法對于解決高維、復雜的RNA生物數據的聚類問題十分有效。
