一種面向SNP選擇的K-Center算法

曹莉敏 周從華

(江蘇大學計算機科學與通信工程學院 江蘇 鎮江 212013)

Information gain

0 引 言

遺傳疾病是由于遺傳物質改變(如基因突變和染色體畸形)而造成的疾病，目前人類已經發現約有4 000種遺傳疾病，以不同程度影響著人類的健康。近些年，隨著DNA微陣列的不斷發展和進步，人類基因層面的研究為找到疾病的發病機理提供了有力的依據。伴隨著人類基因組關聯研究(Genome-Wide Association Study，GWAS)的發展，在精神分類癥、哮喘病等方面的研究得到了良好進展。在全基因組關聯研究中主要利用患病-對照的單核苷酸多態性(Single Nucleotide Polymorphism，SNP)來發現疾病的致病基因，單核苷酸多態性是指DNA序列上核苷酸位點呈現堿基轉換或者顛換，在目前人類基因組的研究范圍內，大約有四百萬個SNP，一般來說，每1 200個堿基中就會有一個SNP。由此可以看出，遺傳疾病的患病與SNP存在及其重要的聯系，并且在巨大而且高維的SNP數據中，使得SNP數據存在著大量的冗余和噪聲，會造成“維度災難”。因此，從SNP數據的特點出發，從大量的SNP數據中盡量地剔除噪聲和冗余，選出具有代表性的SNP子集至關重要。

在進行SNP子集的選擇時,SNP數據存在高維、少樣本的問題，SNP數據之間存在著連鎖不平衡性，使得SNP數據之間存在著關聯性。K-Center算法在聚類時，由于初始聚類中心的確定對于聚類效果來說至關重要，一旦初始聚類中心選擇不好，對聚類結果將產生較大的影響，因此本文利用對稱不確定性來衡量SNP數據之間的相關性，利用初始SNP數據的密度來確定初始聚類中心。(1) K-Center算法在進行聚類時未考慮到數據SNP之間的相關性，將對稱不確定性引入到K-Center中；(2) 在K-Center聚類時，每個SNP數據對于分類準確率貢獻程度的不同，結合SNP數據，采用信息增益來衡量SNP的重要程度，將信息增益引入到密度函數中，來確定高質量的聚類的初始聚類中心。結合上面兩點，提出一種改進的K-Center算法——基于密度的K-MSU算法，利用K-MSU算法將SNP數據聚成K個簇，根據K-MSU算法從K個簇中選出信息SNP子集，將K個子集結合形成最終的信息SNP子集。

1 相關工作

1.1 SNP選擇

SNP的特征選擇問題就是從大量SNP中選出與疾病分類相關的特征子集，然后構建分類模型展開疾病的預測。在目前的特征選擇中有基于機器學習的特征子集選擇方法，大致可以分為過濾法、包裹法以及嵌入式法；如今機器學習和深度學習的發展變得更加火熱，并且給特征選擇提供一個新的思路，能學習到數據中一些有效的表達，得到關鍵的特征[1]。研究者基于以上方法做出很多改進。文獻[2]針對傳統的包裹式方法[3]具有較高的時間復雜度，并且忽略了基因之間的連鎖不平衡性的問題，提出一種基于最大關聯最小冗余(MCMR)的度量方法來構造SNP特征子集，選出具有代表性的SNP子集。文獻[4] 對于提高SNP的準確度和特異性，提出一種基于遺傳算法的SNP特征選擇方法，使得特征數量變少的同時準確率得到提升。文獻[5]指出在高維的全基因組關聯(GWAS)病例對照數據中，存在著較多與疾病無關的單核苷酸多態性(SNP)，提出一種特征子空間的分層抽樣方法，利用隨機森林為SNP數據生成決策樹，獲得更高的準確率和更低的誤差。文獻[6]運用一種新的共享機制，基于搜索策略與聚類方法，利用身高與絕經后乳腺癌和結直腸癌風險之間的聯系，從而可以系統地找到SNP的致病基因，提高SNP預測準確度。通過上述文獻可以看出，選擇高準確率的信息SNP特征的方法應用較為廣泛。

1.2 K-Center算法

K-Center算法是一種對K-Means算法的優化算法，也是一種基于屬性劃分的聚類方法。該算法可以將數據集劃分為簇內相似度高，簇間數據相似度低的簇，但也存在著一些缺點，比如初始聚類中心對聚類結果影響比較大，K值需要事先預定，因此，K-Center和K-Means存在共同問題。文獻[7]提出基于散射中心特征的K-Center算法，將散射中心作為距離計算的特征用于目標的識別，從而提高目標準確識別的精度。文獻[8]提出基于距離系數的K-Center算法，有效避免中心選擇最遠和最近的方法的缺點。文獻[9]則是基于密度的優化初始聚類中心點的選擇的算法，選擇高密度的中心線作為聚類中心，兩種方法都得到了較好的簇集合。文獻[10]在運行速度上提出基于采樣的K-Center MapReduce近似算法，在達到同效果的情況下，提高了運行速度。文獻[11]針對K-Center算法的時間復雜度問題，利用近似線性的逼近算法來找到K-Center中每個簇的大小的界限，改善每個簇的最低和最高界限，減少K-Center算法的時間復雜度。從上述文獻可以看出，K-Center算法對大數據集可以保持其高效性，對序列數據可以有效地進行數據降維，將K-Center算法應用于遺傳疾病的研究相對較少，所以可以將K-Center算法應用于SNP的特征選擇中。

2 基于密度的K-MSU算法的SNP選擇

2.1 K-Center算法原理

K-Means算法屬于劃分層次的聚類方法，從直觀上描述的是簇內樣本的緊密程度，在K-Means聚類中使得同類數據之間的緊密度大，而不同類之間的分離度大。在K-Means的發展中，雖然該算法簡單高效，但也存在著一些問題，比如無法有效地對位點之間的相似進行挖掘，隨機選擇初始聚類中心等方面對聚類結果產生重要影響。K-Center就是K-Means算法的一種改進，區別在于簇中心更新的方式，K-Means算法使用原始的用簇均值作為聚類中心的方法，而K-Center用在這個簇中的某一個樣本點與簇中其他樣本點距離最小的絕對誤差最小的點作為簇中心。

對于給定的數據集D={x1,x2,…,xn}，K-Center算法從數據集D中隨機選擇K個數據點作為初始簇中心點{μ1,μ2,…,μn}；計算樣本xj與簇中心點的距離，根據距離將樣本xj劃入到與之距離最小的一個簇中，更新簇中心；重復上述過程，直到簇中心向量不再更新，聚類結果保持不變，聚類算法停止；得到最終劃分的簇{c1,c2,…,ck}。

K-Center算法的距離計算公式為歐式距離，具體的定義如下：

(1)

式中：n代表總樣本的屬性個數；xi和yi分別表示樣本x和樣本y的第i和第j個屬性。

K-Center算法的目標函數公式，具體定義如下：

(2)

式中：μi為簇中心向量；E為數據樣本到簇中心距離的誤差平方和。

K-Center算法的描述見算法1。

算法1K-Center算法

輸入：數據集D={x1,x2,…,xn}、最大迭代次數n,聚類簇數k。

輸出：聚類得到的k個簇{c1,c2,…,ck}。

1. 從數據集D={x1,x2,…,xn}中隨機選擇k個樣本作為初始簇中心向量：{μ1,μ2,…,μn}

2. Repeat

3. fori=1 tondo

4. 將簇劃分c初始化

5. 根據式(1)計算每個樣本與各簇中心向量的距離

6. 將樣本劃入與之最近的簇中心所在的簇ci

7. end for

8. Forj=1 tondo

9. 利用目標函數計算各個簇，簇內各樣本點距離絕對誤差最小的點|arg minE|作為新的簇中心

10. 將樣本點放入到最近的簇中心所在的簇中

11. end for

until當前的簇中心不再改變，目標函數收斂

2.2 對稱不確定性

信息熵是一種度量樣本集合純度的方法，表示隨機變量的不確定程度，H(X)表示隨機變量的信息熵，給定一個隨機變量X={x1,x2,…,xi}，則信息熵的定義為：

(3)

H(X|Y)代表在隨機變量Y={y1,y2,…,yj}下隨機變量X的條件熵，則條件熵的定義為：

(4)

式中：P(xi)和P(yj)分別代表變量X=xi和Y=yj下的概率；P(xi|yj)代表在隨機變量Y=yj條件下隨機變量X=xi的概率。

信息增益是一種衡量樣本特征重要性的方法，信息增益越大表示該樣本越重要，兩個變量之間的相關性越強，則信息增益的定義如下：

IG(X|Y)=H(X)-H(X|Y)

(5)

在信息增益中，為了使隨機變量之間盡量避免因為變量單位和變量值對結果的影響，采用規范化的信息增益來表示隨機變量X和Y之間的相關性，即對稱不確定性SU(X,Y)，則對稱不確定性的定義如下：

(6)

對稱不確定性的取值區間在(0,1)之間，當SU(X,Y)越大的時候，說明X和Y之間的相關性越大，當SU(X,Y)=0時，說明X和Y之間相互獨立；反之，說明X和Y之間有較強的關聯。通過對稱不確定性可以計算出SNP之間的依賴程度，即兩者之間的相關性。

2.3 K-MSU算法

在K-Center算法中，實現的原理比較簡單，容易實現，算法的可解釋性比較強，其初始聚類中心的選擇，對聚類效果也有較大的影響。在傳統的K-Means算法中并不能挖掘位點之間的相關性，K-Center算法是對K-Means算法的一種改進，但其度量方式也不能挖掘SNP位點之間的相關性?？紤]到SNP位點之間存在連鎖不平衡性，使得SNP數據之間具有強相關性，每個SNP對分類的準確性的重要程度不同，因此首先利用對稱不確定性來確定SNP位點之間的相關性，然后用基于密度和信息增益的方法來選擇SNP中重要的初始聚類中心，結合這兩點提出一種新算法——K-MSU算法。

(1) 距離度量。在SNP數據中，SNP位點之間存在著較強的相關性，因此在原來的K-Center聚類算法距離度量公式中加入對稱不確定性，當SNP位點之間相關性越強時，則SU(X,Y)值越大，此時的歐式距離就會越小，在第一輪的距離度量計算中，每個SNP特征與初始的簇中心的距離度量如下：

(7)

(2) 初始聚類中心的確定。利用密度確定初始聚類中心，在聚類方法中是常用的一種，并且SNP特征對分類結果的重要程度是不同的，因此將密度與信息增益相結合，提出基于信息增益的密度初始聚類中心的選擇方法。設數據集X={x1,x2,…,xi,…，xn}，n為數據的總數，密度值用來衡量初始聚類中心的選擇，對于總的數據樣本中的任意一個樣本數據點xi的密度Dens(xi)，具體定義如下：

(8)

式中：r為給定的有效密度半徑；N為該半徑內所包含的數據點的總數；xj為以xi為圓心r為半徑的第j個數據點；d(xi,xj)為數據點xi和xj改進的歐式距離。

結合SNP數據對最終分類準確率的貢獻程度，利用信息增益增加每個樣本的權值密度來選擇聚類中心點。因此，對于每一個SNP數據點中的樣本點xi的權值wDens(xi)密度，具體的定義如下：

(9)

對于總的數據點來說，其數據點的平均距離定義如下：

(10)

通過計算數據樣本點之間的平均距離，可以有效地反映樣本數據集的整體離散程度，為更好地確定領域半徑提供有效的依據。

鄰域半徑是一種有效密度半徑計算的方法，鄰域半徑不僅參與密度的計算，同時通過鄰域半徑還可以得到在半徑內所包含的數據點數量，合適的鄰域半徑定義如下：

(11)

該鄰域半徑采用的是數據點之間的標準差的方式來確定，在標準差中能夠有效地反映出當前所在的變量和期望變量之間的偏離程度。標準差越小，則鄰域半徑越小，在鄰域半徑內的數據點越緊湊。

根據權值密度確定聚類中心的整體步驟見算法2。

算法2基于信息增益的密度確定聚類中心算法

輸入：數據集S={x1,x2,…,xn}，SNP的特征數n，鄰域半徑r，期望劃分的聚類數目K，最大迭代次數m。

輸出：含有K個對象的初始聚類中心點集T，去除噪聲的數據集P。

1. 令集合D={}，T={}，P={}，q={}。

2. 根據式(7)、式(10)和式(11)計算鄰域半徑r

3. fori=1 tomdo

4. 遍歷集合中的每一個數據對象xi，給xi加入密度標簽，取R=r，根據式(8)計算xi的初始密度值Dens(xi)

5. end for

6. 選取出前2K個初始密度值較大的對象放入到集合D中

7. 通過式(3)-式(5)計算前2K數據樣本的信息增益，存放到信息增益表q中

8. Repeat

9. 取集合D中無中心點或邊界點標簽的第一個對象di

10. fori=1 tomdo

11. 將di作為一個新簇集的聚簇中心點，加入中心點標簽，遍歷集合X中的每一個對象xi，且xi≠di，如果xi是di密度可達點且xi無中心點或邊界標簽，則給xi加入邊界點標簽，納入到該簇內

12. 更新di密度標簽，從q表中取其信息增益值，取R=r，根據式(9)計算di的權值密度wDens(pi)

13. end for

14. until

15. 遍歷集合D中全部對象，即集合D中的每個對象都含有中心點或邊界點標簽

16. 選出集合D中具有中心點標簽且密度值較大的前K個對象加入到集合T中，遍歷集合X中的每一個對象xi，如果xi無中心點或邊界點標簽，則給xi加入噪聲標簽；否則將xi加入到集合P中

(3) 算法步驟。K-MSU算法見算法3。

算法3K-MSU算法

輸入：數據集D={x1,x2,…,xn}、聚類簇數k，最大迭代次數n。

輸出：聚類得到的k個簇{c1,c2,…,ck}。

1. 從數據集D={x1,x2,…,xn}中利用改進的基于信息增益的密度的聚類中心選擇方法選擇初始聚類中心：{μ1,μ2,…,μn}

2. fori=1 tondo

4. 將樣本劃入與之最近的簇中心向量所在的簇ci

5. end for

6. fori=1 tondo

12. 將樣本劃入與之最近的均值向量所在的簇中

13. else

14. 保持當前簇中心不變

15. end for

until當前簇中心不在改變，目標函數收斂

K-MSU利用改進的基于密度的聚類中心選擇方法選擇出初始的聚類中心；利用改進的歐式距離計算每個SNP數據歐式距離，并將其劃分進入相應的簇中，完成第一次迭代；更新聚類中心，利用目標函數計算出簇內各樣本點之間距離絕對誤差最小的點，作為聚類中心；利用當前的聚類中心對簇進行重新劃分，當簇中心不再發生改變，或者迭代次數達到最大，目標函數收斂時，算法完成。

2.4 K-MSU算法在SNP選擇中的應用

本文將K-MSU算法和粒子群算法相結合應用到SNP的特征選擇中。利用K-MSU算法進行特征的聚類，將聚類中的每個簇，利用粒子群優化算法進行特征的選擇，選擇出最佳的SNP特征。

(1) 粒子群算法。粒子群優化算法(PSO)是由Kenned等在1995年提出的一種找尋最優解的算法。該算法是一種群體智能優化算法，主要來源于對鳥類找尋食物的研究，在鳥群搜索的過程中，鳥群都不知道食物的主要位置，每只鳥通過它們之間交流，互通信息，來判斷當前鳥類所找的位置是否是最優解，最終整個鳥群將聚集在食物源的周圍，此時就是整個集群的最優解。PSO中通過位置、速度和適應值三個指標來表示粒子的特征。在粒子的群體運動中，適應度是根據與實際問題相關的適應度函數來進行計算的。以當前選擇的信息SNP來說，最直觀的適應度計算是該信息SNP子集的樣本重構度來作為當前的適應度值。粒子的速度和位置則是通過個體極值和群體極值來更新個體的位置，公式為：

(12)

(13)

式中：vij表示當前粒子的運動速度；wij表示當前粒子的慣性，該參數的取值一般為上一個粒子的運動速度，該參數使運動保持慣性，并且可以增加搜索空間；r1、r2是用來增加粒子多樣性的參數，通過設定一個擾動值，可以防止粒子的飛行超出解的空間范圍；常數c1、c2一般的取值為[0,1]，被稱為學習因子，控制當前粒子靠近最優解和飛行過程的全局最優的速率；pbestij(t)表示在t迭代后所保存的最好的值，此時是一個局部最優值；gbestij(t)表示整個粒子集群的最優值，是一個全局最優解，這兩個值提供了全局最優解的方向；xij表示每個粒子的位置。粒子群算法的流程如下：

① 初始化粒子群，包括樣本的數量，每個粒子的位置和速度；

② 根據適應度函數，計算每個粒子的適應度值；

③ 對每個粒子，用適應度值和個體極值進行比較，當個體極值小于適應度值時，用適應度值替代個體極值；

④ 對每個粒子，用適應度值和全局極值進行比較，當全局極值小于適應度值時，用適應度值替代全局極值；

⑤ 根據公式更新粒子的速度和位置；滿足最大迭代次數則退出，否則返回②。

(2) 基于K-MSU算法的粒子群算法的選擇策略。從聚類后的特征中選擇出合適的特征是進行特征提取和特征選擇的常用方法。一般使用基于過濾式和包裹式的兩種特征選擇方法，但過濾式選擇方法在選擇的時候沒有特定的學習算法進行參與，所以選擇的特征對于目標算法來說并不是最佳的，而基于包裹式的特征選擇方法計算的開銷比較大。本文采用粒子群算法來進行特征選擇。對每個簇進行粒子群算法找尋出最優解，將從K個簇中選擇出來的SNP特征組成新的信息SNP。但是考慮到在不同的簇中，SNP特征的個數也有所不同，如果我們在進行信息SNP特征選擇的時候，將每個簇同等地看待自然是不合理的，因此我們通過簇的大小，來調整每個簇中需要選擇的SNP特征的個數。參數ηk代表當前簇應該選擇的SNP個數，ηk是一個大于等于1的整數，參數ηk的定義如下：

(14)

圖1 SNP特征選擇過程

3 實 驗

3.1 實驗環境與數據

實驗環境：編譯工具Python 3.6.0，操作系統Windows 10，CPU Intel(R)Core(TM)i5-3230M雙核處理器，主頻2.6 GHz，內存8 GB。

實驗數據由無錫市精神衛生中心提供。數據格式為基因型SNP數據，在SNP數據中每個數據都帶有標簽，標明樣本是否患病。SNP數據的具體描述如表1所示。

表1 SNP數據集描述

3.2 評價指標

(1) 聚類效果的評價指標。在目前的聚類評價指標中，輪廓系數結合了聚類的凝聚度和分離度，對聚類效果的評價具有全方位的展現。因此本文中將使用輪廓系數來作為聚類效果的評價指標，通過使用類內距離和類間距離的比值來盡可能評價類內數據之間的緊密度和類間距離的分離度，評價該聚類的效果。類間距離和類內距離，以及輪廓系數的定義如下：

(15)

(16)

(17)

式中：類內距離a(i)是樣本xi與同類中其他樣本的平均值，ni代表樣本xi所在類別的總數目；類間距離bi為樣本xi到其他每個類中樣本平均距離的最小值，pi為樣本xi所在的類別，mp為該類別的第p類數據的聚類中心，K為劃分簇數；S(i)為劃分為k類的聚合效果指標，S(i)越大說明聚類效果越好。

(2) 分類情況的評價指標。評價標準使用精準率(precision)、召回率(recall)、F1值和預測準確率(Acc)來進行評價，其最終的評價標準可以由“混淆矩陣”進行計算，如表2所示。

表2 “混淆矩陣”結果

通過“混淆矩陣”，可以得出評價分類效果的指標。

(18)

(19)

(20)

(21)

3.3 數據預處理

由于原始的數據類型是SNP的基因型格式，在后續的實驗中不利于聚類的實驗，因此在實驗之前需要將SNP的基因型格式編碼，將基因型AA、Aa、aa分別用0、1、2表示。編碼完成后，對缺失值使用K近鄰進行填充。由于編碼后的SNP數據會存在有大量缺失情況，因此將大于10%的缺失值樣本剔除。

3.4 實驗結果與分析

(1) 信息SNP選擇實驗及分析。信息SNP選擇主要采用三種方法對比，分別是K-Means、K-Center和K-MSU，對比采用的評價指標為利用聚類后使用粒子群算法選擇出的信息SNP子集對結果的平均預測準確率。預測準確率可以作為信息SNP子集是否可以代表其余所有SNP的重要信息度量，預測準確率越高，說明信息SNP越可以代表剩下的非信息SNP。在這兩個數據集上，用三種方法對這兩個數據集進行聚類實驗，對K個簇分別進行粒子群算法選擇信息SNP，最后采用KNN對信息SNP預測，得到預測準確率，每組實驗重復三次對結果取平均值，實驗結果如圖2和圖3所示。

圖2 在數據集Dataset1上的平均預測準確率

圖3 在數據集Dataset2上的平均預測準確率

可以看出，算法預測準確率隨著聚類簇個數K的增加而增大，當簇的個數為8時，平均預測準確率最高。利用K-MSU進行聚類，由粒子群算法進行特征選擇的信息SNP達到的平均預測準確率最高。在接下來的實驗中，將采用簇個數為8時選擇的信息SNP數據進行分類實驗。

(2) 聚類效果評估。本次實驗采用輪廓系數來進行聚類效果的評估，當輪廓系數越接近于1時類內距離越小，類間距離越大，說明聚類效果越好。三種算法在數據集Dataset1和Dataset2上的聚類評價結果如表3和表4所示。

表3 數據集Dataset1下的三種算法的聚類評價結果

表4 數據集Dataset2下的三種算法的聚類評價結果

可以看出，K-MSU算法在兩個數據集上有較好的聚類效果，在7組實驗中分別都有5次和6次達到最好。對K-MSU算法而言，當聚類個數在8個時，兩個數據集的效果最好，后續將使用K=8作為最終的聚類數目。

(3) SNP子集分類實驗及分析。在信息SNP子集選擇后，為了進一步驗證信息SNP選擇的重要性，采用K-MSU算法、K-Center算法[6]、K-Means/粒子群算法[9]、K-Center/粒子群算法[12]、K-MSU/粒子群算法、ReliefF[13]和MCMR[2]算法對信息SNP進行篩選，將選擇出來的信息SNP子集利用SVM、DT和神經網絡作為分類模型，進行分類實驗。此次分類實驗的主要評價指標為F1值和分類準確率Acc。實驗結果如表5所示。

表5 利用不同分類器進行SNP子集評價的結果

對比K-MSU算法和K-Center算法可以得出，改進的K-MSU算法的F1值和分類準確率都有所提高；在不同的分類實驗中，K-MSU/粒子群算法與K-Means粒子群算法、K-Center/粒子群算法進行篩選出的SNP子集在不同的分類器上表現分類效果是不同的，與兩個常用的特征選擇方法ReliefF和MCMR方法相比具有更好的分類結果。

4 結 語

SNP數據具有高維少樣本，SNP數據之間存在強相關性的問題，且K-Center方法在初始聚類中心的選擇中沒有好的方法，結合這兩個問題,本文提出一種基于K-Center方法改進的K-MSU方法并將其應用到SNP的特征選擇中。首先對SNP數據進行初步的篩選，再利用K-MSU方法對SNP數據進行聚類，最后對聚類后的每個簇使用粒子群算法構造信息SNP子集。在最后的分類結果中，對于SNP子集的特征選擇均表現出有效性。本文的亮點：(1) 對于聚類中心的確定時，不采用隨機選取，而是利用信息增益和密度選擇重要程度相對高的SNP數據點作為聚類中心；(2) 在對SNP聚類時，由于原始的K-Center方法無法挖掘出位點之間的關聯性問題，因此引入了對稱不確定性，來衡量位點之間的關聯性，將其應用到K-Center中。未來將對其分類方法進行改進，提高SNP的分類準確率。