邢 斌 周從華 張付全 張 婷 蔣躍明

（1.江蘇大學(xué)計(jì)算機(jī)科學(xué)與通信工程學(xué)院 鎮(zhèn)江 212013）（2.無錫市精神衛(wèi)生中心 無錫 214151）（3.無錫市婦幼保健院 無錫 214002）（4.無錫市第五人民醫(yī)院 無錫 214073）

1 引言

遺傳病是由致病基因所控制的疾病，這種由于遺傳物質(zhì)發(fā)生改變引起的疾病病種多，目前已發(fā)現(xiàn)的遺傳病超過3000種，而且具有發(fā)病率高、先天性、終生性、家族性等特點(diǎn)，對(duì)人類健康產(chǎn)生巨大影響。近年來隨著世界范圍內(nèi)人類基因組研究（Ge?nome-Wide Association Study，GWAS）開展實(shí)施，基因測(cè)序技術(shù)依靠對(duì)遺傳特征的有效挖掘，在疾病診斷、預(yù)測(cè)和治療等方面發(fā)揮著更加重要作用［2］，而GWAS是在單核苷酸多態(tài)（Single Nucleotide Poly?morphism，SNP）性的基礎(chǔ)上展開研究的，單核苷酸多態(tài)性主要是指基因組水平上的單核苷酸變異導(dǎo)致脫氧核糖核酸（DNA）序列多態(tài)性的現(xiàn)象，是人類最常見的一種遺傳變異，研究表明SNP對(duì)人類疾病有直接或間接的聯(lián)系，并且SNP具有穩(wěn)定性好、頻率高、采集容易等優(yōu)點(diǎn)，因此對(duì)SNP數(shù)據(jù)進(jìn)行研究有重要意義。在此背景下，許多機(jī)器學(xué)習(xí)算法在SNP數(shù)據(jù)中得到了廣泛的應(yīng)用。但是由于SNP數(shù)量過大，維數(shù)過高，數(shù)據(jù)中存在冗余和噪聲，研究中必須要考慮“維數(shù)災(zāi)難”問題。因此，SNP數(shù)據(jù)分析的初始階段一般是選擇SNP中信息量最大的子集即信息SNP子集，以提高算法的性能，降低時(shí)間要求。

信息SNP子集的選擇本質(zhì)上是特征選擇，特征選擇（Feature Selection，F(xiàn)S）是在保持原始數(shù)據(jù)準(zhǔn)確表示的同時(shí)，顯著降低特征空間的維數(shù)的過程。然而由于SNP數(shù)據(jù)之間的非相互獨(dú)立的特點(diǎn)，現(xiàn)有的特征選擇方法難以挖掘SNP位點(diǎn)之間的相關(guān)性，從而漏掉重要的遺傳信息，最終降低算法的分類效果。鑒于上述問題，本文采用基于信息論的新的相似度度量方法，并將其應(yīng)用到K-means中，把SNP依據(jù)一定的關(guān)系劃分為若干簇，然后使用粒子群算法從每個(gè)簇中選出一個(gè)或多個(gè)能夠代表整個(gè)簇的信息SNP，構(gòu)造出最終的信息SNP子集。

2 相關(guān)工作

2.1 SNP選擇研究現(xiàn)狀

目前SNP選擇主要有兩種方法，一種是基于統(tǒng)計(jì)檢驗(yàn)的關(guān)聯(lián)研究，分兩個(gè)步驟，第一步，通過生物實(shí)驗(yàn)技術(shù)在全基因組上掃描篩選有效位點(diǎn)，第二步查驗(yàn)SNP位點(diǎn)基因分型，通過關(guān)聯(lián)分析識(shí)別信息SNP；另一種是基于機(jī)器學(xué)習(xí)的SNP選擇，這種SNP選擇本質(zhì)上是特征選擇問題，可分為過濾式和封裝式兩種。過濾式設(shè)置評(píng)價(jià)指標(biāo)給每個(gè)SNP打分，優(yōu)點(diǎn)是計(jì)算量小，但忽略了SNP之間的內(nèi)在聯(lián)系，不容易得到最優(yōu)的特征子集；而包裹式方法雖然計(jì)算量大，卻可以把學(xué)習(xí)器和評(píng)價(jià)指標(biāo)結(jié)合起來，最終得到最優(yōu)解。近年來，許多人嘗試對(duì)這兩種方法進(jìn)行改進(jìn)，Raid Alzubi等［1］將過濾式和包裹式方法結(jié)合在一起提出一種混合特征選擇方法來檢測(cè)信息SNP并選擇最優(yōu)SNP子集，取得了不錯(cuò)的效果。為解決SNP數(shù)據(jù)高維問題，Cong等［2］提出了一種基于主成分分析的基因數(shù)據(jù)降維算法，用于低維空間SNP位點(diǎn)的聚類。上述兩篇文獻(xiàn)雖然在結(jié)果上較前人有一定的進(jìn)步，但是仍忽略了SNP與SNP子集之間的相關(guān)性問題，Liao［3］將信息論引入到SNP選擇中，用于衡量SNP子集之間的相關(guān)性，該方法顯著提高了標(biāo)簽SNP選擇的效率和預(yù)測(cè)精度，但該方法沒有將SNP之間的相似度和聚類有效結(jié)合在一起。

2.2 K-means聚類算法

K-means算法是一種迭代求解的無監(jiān)督聚類算法，通過選取K個(gè)對(duì)象作為初始聚類中心計(jì)算每個(gè)樣本點(diǎn)與聚類中心的距離即相似度，將每個(gè)樣本點(diǎn)分配到最近的聚類中心，最終使得相同簇內(nèi)的樣本點(diǎn)相似度盡可能高，不同簇間的相似度盡可能低。

2.2.1 歐氏距離

歐氏距離是衡量兩個(gè)同一樣本集樣本對(duì)象差異性，距離越大樣本差異度越大。歐氏距離公式如下：

式中，xiμ為樣本對(duì)象xi的第μ個(gè)屬性。

2.2.2 簇內(nèi)誤差平方和

在歐氏距離基礎(chǔ)上，進(jìn)行變形，將其中一個(gè)樣本對(duì)象換為簇中心得到簇內(nèi)誤差平方和，簇內(nèi)誤差平方和用來度量聚類效果的好壞。

3 基于改進(jìn)的K-means算法的SNP選擇方法

3.1 傳統(tǒng)K-means在SNP選擇上的缺陷

傳統(tǒng)的K-means聚類使用SNP之間的歐氏距離來度量SNP之間的相似度，這并不能挖掘出SNP位點(diǎn)之間生物學(xué)上的聯(lián)系性，鑒于上述問題［1］將互信息引入到聚類中，用作SNP相似度度量，雖然起到很好的效果，但仍沒解決傳統(tǒng)K-means算法中另一缺陷，只考慮兩個(gè)SNP位點(diǎn)之間的相似度，然而在實(shí)際SNP數(shù)據(jù)集中，單個(gè)SNP往往和某個(gè)SNP子集有著強(qiáng)關(guān)聯(lián)性，鑒于上述兩個(gè)問題，提出一種新的基于信息論的SNP選擇算法K-MIGS。

3.2 基于信息論的相似度度量

3.2.1 信息熵與互信息

信息熵用來衡量信息的不確定性，單個(gè)SNP位點(diǎn)的信息熵可表示為

假設(shè)樣本S中某一SNPX有T個(gè)可能取值，其中p(xi)表示位點(diǎn)X中第i類樣本出現(xiàn)的頻率。Tx表示位點(diǎn)X所有屬性的個(gè)數(shù)。聯(lián)合信息熵：

則SNP位點(diǎn)X和Y的互信息可表示為

我們使用兩個(gè)位點(diǎn)的互信息來衡量它們之間的相似度，互信息越大，表示相似度越大，每個(gè)特征與初始簇中心的距離度量公式表示如下：

接下來我們?cè)噲D利用MI(S1;SNPX)和MI(S2;SNPX)來構(gòu)建SNPx與SNP子集之間的相似度，有下式：

假設(shè)n個(gè)SNP位點(diǎn)上存在m種，那么SNP子集內(nèi)聯(lián)合熵可表示為

但是如果S1、S2兩個(gè)子集SNP個(gè)數(shù)是不相等的，那么I(S1;SNPX)和I(S2;SNPX)將沒有可比性，因?yàn)閱蝹€(gè)SNP和大的SNP子集會(huì)得到大的互信息值，為了解決這個(gè)問題，我們結(jié)合了互信息和信息熵提出了一種單個(gè)SNP和SNP子集的相似度度量方式MIGS，公式如下：

其中，H(S)、H(SNPX)分別表示SNP子集S和單個(gè)SNP的信息熵，MI(S;SNPX)表示它們之間的互信息，考慮到SNP位點(diǎn)之間存在強(qiáng)相關(guān)性的特性，我們?cè)谠璌-means聚類的距離度量歐氏距離中引入互信息的概念，則第一輪迭代計(jì)算中，每個(gè)特征與初始簇中心的距離度量公式表示如下：

其中MId(xi，uj)表示特征xi與初始簇中心uj的相似度，表示特征與初始簇中心的歐氏距離，MIGS(xi，Ci)表示特征x i與初始簇i之間的相似度。

3.2.2 簇中心的更新

在K-MIGS中，初始的簇中心的選擇與傳統(tǒng)K-means一致，不同的是在簇中心更新時(shí)，傳統(tǒng)K-means采用均值向量的方式更新簇中心，這并不適合用在K-MIGS算法中，因此我們采用新的簇更新方式，選取與均值向量最近的一個(gè)SNP作為簇中心。

3.3 算法K-MIGS步驟

結(jié)合章節(jié)3.1和3.2，則算法K-MIGS的整體步驟如算法1所示。

算法1：K-MIGS算法

輸入：數(shù)據(jù)集D={x1，x2，…，xm}聚類簇?cái)?shù)kMAX_N

1）從數(shù)據(jù)集D={x1，x2，…，xm}中隨機(jī)選擇k個(gè)樣本作為初始均值向量(μ1，μ2，…，μk)

2）for j=1 to k do

3） forj=1 j=1 to mdo

4）若簇CiSNP數(shù)為1，根據(jù)式（6）、（7）計(jì)算與μj的相似度距離，否則根據(jù)式（10）、（11）計(jì)算SNP子集S與SNPi的相似度距離MId(S1，SNPi)

5） end for

6）end for

7）for i=1 to m do

8）將xi劃入與其相似度距離最小的簇

9）end for

10）repeat

11）for i=1 to k do

12）計(jì)算新的均值向量，并通過式（1）計(jì)算x( )x∈Ci與μ'i的距離

13）選擇最小的d(x，μ'i)作為簇Ci新的簇中心

14）end for

15）until迭代次數(shù)達(dá)到閾值或簇中心不再更新

其中，2～6行根據(jù)式（6）、（9）、（10）分兩種情況計(jì)算SNP之間的相似度距離，分為當(dāng)前簇SNP數(shù)為一和SNP數(shù)不為一；7～9行將所選SNP劃分到相似度距離最小的簇，11～13行根據(jù)式（1）將簇中心更新為與均值向量歐氏距離最小的SNP。

3.4 K-MIGS算法在SNP選擇中的應(yīng)用

本節(jié)將結(jié)合K-MIGS算法和粒子群算法對(duì)每個(gè)簇中的SNP進(jìn)行選擇。

粒 子 群 算 法（Particle Swarm Optimization，PSO）是一種群體智能算法，該算法模仿鳥群群體合作來最快獲得食物的基本原理。在粒子群算法中，每一個(gè)候選解都可以看作是“鳥群中的個(gè)體”，即一個(gè)粒子。在粒子運(yùn)動(dòng)過程中每個(gè)粒子根據(jù)自身經(jīng)驗(yàn)和相鄰粒子經(jīng)驗(yàn)，通過對(duì)速度進(jìn)行改變來調(diào)整自己的位置。粒子通過最佳粒子流在問題空間內(nèi)運(yùn)動(dòng)，迭代一定次數(shù)或者達(dá)到給定的最小誤差得到最優(yōu)解。

將粒子群算法應(yīng)用到SNP選擇中包括兩個(gè)主要部分，種群初始化和粒子更新。

我們?yōu)镾NP選擇問題設(shè)計(jì)了種群：

其中n代表種群P的大小，m代表SNP的數(shù)量，pij=1表示選擇了第i個(gè)粒子的第j個(gè)SNP。

通常，在種群的初始化過程中每個(gè)SNP的選擇為信息SNP的概率由確定，其中k代表選擇的SNP的數(shù)量，m代表SNP的總數(shù)。在每次迭代中依據(jù)pbest和gbest的參數(shù)對(duì)粒子進(jìn)行更新，每個(gè)粒子根據(jù)以下方程式更新：

其中w是慣性權(quán)重，a1和a2是加速學(xué)習(xí)因子，r1和r2為0～1的隨機(jī)數(shù)，和為更新前和更新后的粒子i速度矢量的第j維分量，為更新前粒子i位置矢量的第j維分量，pbestij=(pi1，pi2，…，pij)表示粒子i個(gè)體經(jīng)歷過的最好位置，gbestj=(g1，g2，…，gj)表示種群所經(jīng)歷過的最好位置。

將實(shí)驗(yàn)所需原始數(shù)據(jù)進(jìn)行預(yù)處理，并使用K-MIGS算法進(jìn)行聚類，最后使用粒子群算法從每個(gè)簇中選擇得到最終的信息SNP。

4 實(shí)驗(yàn)

4.1 實(shí)驗(yàn)環(huán)境和數(shù)據(jù)

實(shí)驗(yàn)環(huán)境：編譯工具Python 3.7.3，操作系統(tǒng)Windows10 64位，處 理 器Intel（R）Core（TM）i7-6700HQ，CPU@2.60GHz 2.59GHz，GPU Nvidia GeForce GTX 1060 3G，運(yùn)行內(nèi)存16G，硬盤1.25T。

實(shí)驗(yàn)數(shù)據(jù)集：實(shí)驗(yàn)數(shù)據(jù)來自于無錫精神衛(wèi)生中心，包括兩種SNP數(shù)據(jù)集EN1000（9445SNPS）、E144（2514825SNPS），另外還包含樣本的基因型，患病與否的標(biāo)記，具體描述如表1所示。

表1 數(shù)據(jù)集描述

4.2 實(shí)驗(yàn)評(píng)價(jià)指標(biāo)

1）SNP選擇評(píng)價(jià)指標(biāo)

本文使用SNP選擇常用評(píng)價(jià)指標(biāo)，信息SNP子集對(duì)非信息SNP子集重構(gòu)準(zhǔn)確度來評(píng)價(jià)SNP選擇效果的好壞。

其中pi為信息SNP子集對(duì)非信息SNP位點(diǎn)i的重構(gòu)準(zhǔn)確度，||O是非信息SNP的數(shù)量。重構(gòu)準(zhǔn)確度越高，SNP選擇效果越好。

2）分類預(yù)測(cè)評(píng)價(jià)指標(biāo)

本文使用分類實(shí)驗(yàn)常用的評(píng)價(jià)指標(biāo)預(yù)測(cè)的準(zhǔn)確率（Accuracy）及F-Measure對(duì)分類結(jié)果進(jìn)行評(píng)價(jià)，分類結(jié)果的“混淆矩陣”如表2所示。

表2 預(yù)測(cè)類別與實(shí)際類別的“混淆矩陣”

根據(jù)上表分類結(jié)果評(píng)價(jià)指標(biāo)可由式（14～17）計(jì)算：

4.3 數(shù)據(jù)預(yù)處理

SNP數(shù)據(jù)預(yù)處理分為兩個(gè)子階段，數(shù)據(jù)編碼和數(shù)據(jù)更新，原始SNP具有三種基因型，純合子基因型（AA），雜合子基因型（Aa）以及純合變異基因型（aa），這種基因型數(shù)據(jù)不利于后續(xù)聚類操作，需要對(duì)數(shù)據(jù)進(jìn)行編碼，分別將AA、Aa、aa編碼為0、1、2；數(shù)據(jù)編碼后對(duì)數(shù)據(jù)進(jìn)行更新，更新又包括缺失數(shù)據(jù)填充及不符合標(biāo)準(zhǔn)數(shù)據(jù)的刪除。依據(jù)上述原則處理后，數(shù)據(jù)集G1000和E144分別剩余9298和214935條有效的SNP。

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

4.4.1 聚類實(shí)驗(yàn)及分析

此部分實(shí)驗(yàn)包括四種算法比較實(shí)驗(yàn)，分別為K-means、特征加權(quán)K-means、模糊聚類算法FCM和K-MIGS，主要評(píng)價(jià)指標(biāo)為最后得到不同簇下信息SNP子集的重構(gòu)準(zhǔn)確度。在數(shù)據(jù)集G1000和E144上，分別用四種算法進(jìn)行聚類實(shí)驗(yàn)，使用粒子群算法對(duì)每個(gè)簇進(jìn)行信息SNPs提取，最后計(jì)算信息SNP子集對(duì)非信息SNP子集的重構(gòu)度。實(shí)驗(yàn)結(jié)果如圖1和圖2所示。

圖1 G1000上算法選出的信息SNP對(duì)非信息SNP的重構(gòu)度圖

圖2 E144上算法選出的信息SNP對(duì)非信息SNP的重構(gòu)度圖

由圖1、2可看出，使用K-MIGS/粒子群算法最終提取的信息SNP在兩個(gè)數(shù)據(jù)集上對(duì)非信息SNP具有更好的重構(gòu)度，并且在聚類簇?cái)?shù)為8時(shí)重構(gòu)度達(dá)到最大值，因此后續(xù)分類實(shí)驗(yàn)采用簇?cái)?shù)為8。

4.4.2 分類實(shí)驗(yàn)及分析

分類實(shí)驗(yàn)的目的是更進(jìn)一步檢測(cè)K-MIGS/粒子群算法所選擇的信息SNP子集包含信息的重要程度，此實(shí)驗(yàn)中，采用K-means/粒子群算法、K-MIGS/粒子群算法、特征加權(quán)K-means/粒子群算法（FW-K-means/粒子群）、ReliefF和MCMR算法進(jìn)行信息SNP子集的篩選，使用SVM、DT和神經(jīng)網(wǎng)絡(luò)（Neural Networks，NN）為分類器，主要評(píng)價(jià)指標(biāo)為分類準(zhǔn)確率Acc和F1-measure。實(shí)驗(yàn)結(jié)果如表3所示。

表3 不同的分類器進(jìn)行SNP子集評(píng)價(jià)的結(jié)果

5 結(jié)語

針對(duì)SNP本身具有的高維少樣本特性以及自身存在的遺傳規(guī)律，本文提出了一種基于K-means的算法K-MIGS并將其應(yīng)用到SNP選擇中。對(duì)數(shù)據(jù)進(jìn)行預(yù)處理后，使用提出的方法對(duì)SNP數(shù)據(jù)進(jìn)行聚類，最后使用粒子群算法構(gòu)造最終的SNP子集。在聚類效果評(píng)估和分類實(shí)驗(yàn)中均表明，該方法很大地提升了SNP選擇的有效性。本文提出的方法相比其他SNP選擇方法優(yōu)勢(shì)在于同時(shí)考慮了單個(gè)SNP與SNP子集的相似度對(duì)聚類結(jié)果的影響。本文的后續(xù)工作主要有兩點(diǎn)：一是繼續(xù)優(yōu)化算法以減少單個(gè)SNP與SNP子集的相互關(guān)聯(lián)時(shí)引入的額外時(shí)間復(fù)雜度；二是繼續(xù)優(yōu)化粒子群算法，使其在每個(gè)簇中選擇出信息量更大的SNP。