顧 鑫 陳偉鶴 張付全 張 婷 蔣躍明

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

1 引言

遺傳疾病是指由于遺傳物質(zhì)發(fā)生改變而引發(fā)的疾病，目前遺傳疾病主要分為復(fù)雜疾病和單基因疾病兩種，復(fù)雜疾病主要包括精神分裂癥和哮喘病等，主要是因為基因中的多個單核膽酸多態(tài)性之間的相互作用而形成的，而單基因疾病則是遵循孟德爾遺傳定律。

近年來，隨著DNA 微陣列技術(shù)的不斷進(jìn)步，作為檢測人類常見疾病的遺傳變異的工具，基因組范圍的關(guān)聯(lián)研究也受到了很大的關(guān)注。單核苷酸多態(tài)性（Single Nucleotide Polymorphism，SNP）是指基因組水平上由于單核苷酸變異所引起的序列多態(tài)性。SNP 擁有數(shù)量多、分布范圍廣和穩(wěn)定度高等特點(diǎn)，常被用于復(fù)雜性狀的疾病、群體的基因識別和遺傳解剖等方面的研究，因此，SNP 已經(jīng)成為第三代遺傳標(biāo)記。對SNP的廣泛研究，使得像類風(fēng)濕關(guān)節(jié)病和精神分裂癥等疾病的研究取得了良好的進(jìn)展［1］。然而，大量研究發(fā)現(xiàn)兩個無關(guān)個體的99.9%的基因組序列是相似的。而剩余的0.1%的差異是導(dǎo)致人體產(chǎn)生疾病的關(guān)鍵所在［2～3］，因此對冗余的SNP 進(jìn)行篩選，即從大量的SNP 中選擇具有代表性的信息SNP成為一個重要的課題。

不同遺傳標(biāo)記之間存在非隨機(jī)組合的現(xiàn)象，例如多代遺傳中的SNP，即標(biāo)記不是完全獨(dú)立的。這種現(xiàn)象通常存在于各種物種中，我們將這種現(xiàn)象稱為連鎖不平衡（LD）。在對SNP 進(jìn)行篩選時，考慮到SNP維度較高以及SNP之間存在連鎖不平衡性的特點(diǎn)，使得傳統(tǒng)的機(jī)器學(xué)習(xí)方法在解決它時難免會遺漏掉許多內(nèi)在的遺傳信息。針對上述場景，結(jié)合上述SNP 的特點(diǎn)，本文提出一種基于改進(jìn)蟻群算法的SNP選擇方法，設(shè)計合理的路徑選擇函數(shù)和信息素更新機(jī)制，同時將連鎖不平衡性引入蟻群算法。

2 相關(guān)工作

目前，國內(nèi)外的相關(guān)研究主要是通過生物實驗的方法從樣本中獲取SNP的原始數(shù)據(jù)，如果單純采用生物學(xué)的方法進(jìn)行基因分型，將會面臨消耗時間長，代價昂貴以及難以滿足生物分析數(shù)據(jù)的要求［4～5］。當(dāng)前的SNP 選擇方法有很多，比較成熟的有兩類：基于單塊的方法和基于單體型重構(gòu)的信息SNP方法。

2.1 基于單體型塊的方法

考慮到理論數(shù)量遠(yuǎn)大于人類單體型數(shù)量的基本事實，通過設(shè)定一個評價指標(biāo)來衡量每個SNP，將基因組序列數(shù)據(jù)分成多個離散的單元型塊，然后根據(jù)相應(yīng)的規(guī)則在每個塊中選擇相應(yīng)的信息SNP。Patil 首先提出了使用貪心算法來劃分奇異塊的想法［6］。Chang 等提出了混合貪婪-劃分樹的方法，該方法引入了分支算法定界的思想，一個原始信息SNP選擇問題被劃分為多個獨(dú)立的子問題，最后構(gòu)建出貪婪劃分樹［7］。Liao 提出一種多次蟻群算法選擇SNP集合，通過計算復(fù)雜度和噪聲影響同時提高劃分準(zhǔn)確率，試驗結(jié)果表明該方法有一定改進(jìn)［8］。Prathibh 提出了一種基于遺傳算法的特征選擇算法，該算法減少了特征數(shù)量，提高了基因/SNP組的特異性［9］。

2.2 基于單體型重構(gòu)的方法

Bafna［10］及Halldorsson［11］首先提出了一種基于單倍型重構(gòu)的SNP位點(diǎn)選擇方法。Halperin［12］等描述了一種用于SNP 預(yù)測和選擇的新方法STAMPA，該方法可以不必提前執(zhí)行區(qū)塊劃分，因此使得方法的應(yīng)用范圍更加廣闊。Lee［13］提出了基于SNP 之間的條件獨(dú)立性來標(biāo)記SNP選擇的方法，通過構(gòu)建貝葉斯網(wǎng)絡(luò)，選擇獨(dú)立和預(yù)測性高的SNP的一個子集。Ilhan［14］提出采用克隆選擇算法選擇SNP 子集，其中SNPs 之間的相似性關(guān)聯(lián)被用作其余SNP的預(yù)測方法，能夠更快的識別SNP。Alzubi［15］提出了基于條件互信息最大化和支持向量機(jī)特征遞歸消除融合的混合特征選擇方法，取得了較高的重構(gòu)準(zhǔn)確度。

3 基于蟻群算法的SNP選擇方法

根據(jù)SNP 數(shù)據(jù)分布特點(diǎn)和信息SNP 選擇的難點(diǎn)，本文提出基于改進(jìn)蟻群算法的信息SNP 選擇方法，設(shè)計合理的路徑選擇函數(shù)和信息素更新機(jī)制。為了避免信息素的過分累積，從而引發(fā)局部最優(yōu)，提出對信息素進(jìn)行揮發(fā)，同時將連鎖不平衡性引入蟻群算法的啟發(fā)式函數(shù)，從而對SNP選擇方法進(jìn)行改進(jìn)。

3.1 基于改進(jìn)蟻群算法的SNP子集構(gòu)造

3.1.1 連鎖不平衡

目前已有多種連鎖不平衡性的測量指標(biāo)，包括兩位點(diǎn)和多位點(diǎn)。以下為二等位基因位點(diǎn)的連鎖不平衡度量方法，假定兩個SNP位點(diǎn)的四種數(shù)據(jù)頻率分別為f11，f12，f21，f22，它們滿足式（1）：

其中，f1+標(biāo)識第一個等位基因為1 的單體型的頻率之和，那么連鎖不平衡度量方法如式（2）所示：

兩點(diǎn)連鎖度量D 值范圍太大［16］，具有相似連鎖分布的位點(diǎn)組合之間的D 值將會變得更大。改進(jìn)的策略是對該值作歸一化后再進(jìn)行度量，如式（3）所示：

此時r2的取值位于0～1 之間。研究發(fā)現(xiàn)，r2能直接反應(yīng)關(guān)聯(lián)研究效能，因此當(dāng)前大多方法采用它作為選擇代表性SNP的依據(jù)。

3.1.2 蟻群算法基本原理

蟻群算法最初由意大利學(xué)者Dorigo M。于1991 年首次提出，它本質(zhì)上是一個由仿生學(xué)計算構(gòu)建的群優(yōu)化系統(tǒng)。蟻群算法具有天然的分布式計算機(jī)制、較強(qiáng)的魯棒性和易與其他優(yōu)化算法結(jié)合的特點(diǎn)。

在覓食的過程中，螞蟻將會根據(jù)信息素的濃度來決定移動的方向。因此，當(dāng)環(huán)境中沒有信息素時，螞蟻的行為將是完全隨機(jī)的。而在接下來的過程中，一條路徑上經(jīng)過的螞蟻越多，那么這條路徑上積累的信息素也就越多，之后的其他螞蟻因此更有可能選擇這條路徑，該過程逐漸由隨機(jī)行為轉(zhuǎn)變?yōu)橹悄苄袨椤?/p>

3.1.3 路徑選擇函數(shù)

信息SNP 選擇問題中候選子集的質(zhì)量取決于兩個因素：信息SNP 數(shù)量和信息SNP 對非信息SNP 重新構(gòu)造的準(zhǔn)確度。本文將這兩種因素放在路徑選擇及信息素更新過程中。螞蟻的路徑選擇采用概率機(jī)制，當(dāng)前人工螞蟻選擇下一節(jié)點(diǎn)的概率如式（4）所示：

其中τi和ηi分別表示節(jié)點(diǎn)信息素濃度和節(jié)點(diǎn)當(dāng)前的啟發(fā)信息，而α和β分別為信息素和啟發(fā)因子的權(quán)重參數(shù)，R 表示本次迭代過程中沒有被選中的SNP 位點(diǎn)?？梢酝ㄟ^調(diào)整權(quán)重的方式更新選擇機(jī)制，如當(dāng)α＞β時，表示蟻群選擇路徑更側(cè)重于參考SNP位點(diǎn)上的信息素濃度。

3.1.4 信息素更新機(jī)制

信息SNP 選擇問題中信息SNP 數(shù)量大小可以類比于傳統(tǒng)螞蟻在覓食過程中走過的路程，路程的長度越短則表明SNP 該條路徑越優(yōu)秀。在相同的重新構(gòu)造準(zhǔn)確度下，信息SNP 的數(shù)量越小，SNP 子集越好，位點(diǎn)上的信息素類似于自然界中螞蟻留下的化學(xué)物質(zhì)。信息素累積函數(shù)如式（5）所示：

其中，Tk(t)表示第k只螞蟻在第t 次迭代過程中所構(gòu)造的候選信息SNP 子集，Q為原始數(shù)據(jù)集中的所有SNP 數(shù)量。Lk(t)表示路徑的長度，即子集中包含的SNP數(shù)量。

信息SNP 選擇問題也與天然螞蟻覓食過程不同，信息SNP選擇過程是跳躍的，相反，螞蟻可以從當(dāng)前節(jié)點(diǎn)跳轉(zhuǎn)到任何其他節(jié)點(diǎn)，并且構(gòu)造的候選子集具有無序性，即選擇SNP的順序不影響子集的質(zhì)量。

為了防止信息素在某些位置疊加并引起局部最優(yōu)，必須適當(dāng)削弱位點(diǎn)的信息素。本文模擬自然環(huán)境中的空氣流動，引入了信息素?fù)]發(fā)機(jī)制，通過式（6）實施螞蟻和信息素蒸發(fā)添加新信息素：

?τi( )

t表示迭代后式（2）中在所具有的信息素，并且被設(shè)置為路徑上的總信息素的初始值。

信息素?fù)]發(fā)系數(shù)ρ是螞蟻留在路徑上的信息素的持久權(quán)重，揮發(fā)系數(shù)越小表示信息素每次迭代過程中損失越小，信息素越不容易消失，留下的信息素越多。通過調(diào)整揮發(fā)系數(shù)，可以在一定程度上降低信息素的過度累積，從而盡量避免蟻群算法在后期陷入局部最優(yōu)?？紤]到該算法在初始階段較差的尋優(yōu)能力，在得到更好的解后引起局部最優(yōu)，因此提出一種新的信息素?fù)]發(fā)因子。具體調(diào)整如式（7）所示：

式（7）中，NC_MAX 表示蟻群算法的最大迭代次數(shù)，為固定值。NC 表示蟻群算法在本輪的迭代次數(shù)。在迭代初期，通過設(shè)定較小的信息素?fù)]發(fā)系數(shù)可以加快收斂速度，而在算法的中期和末期，適當(dāng)增加信息素?fù)]發(fā)系數(shù)的值避免局部最優(yōu)。

3.1.5 啟發(fā)式函數(shù)

信息SNP 子集的優(yōu)劣可以由重構(gòu)準(zhǔn)確度來測量，其中Ci為兩位點(diǎn)連鎖不平衡度量，如式（8）所示：

本文中兩位點(diǎn)連鎖測量使用式（3）的r2。其中表示兩個位點(diǎn)i 和j 之間的連鎖值，而n表示已經(jīng)選中位點(diǎn)的數(shù)目，Ci為位點(diǎn)i和已經(jīng)選中位點(diǎn)的連鎖值的平均值，介于0～1 之間。值越大，該位點(diǎn)具有更強(qiáng)的連鎖性，則更有利于成為信息SNP。

Pi表示信息位點(diǎn)i 對其它所有非信息SNP 位點(diǎn)重構(gòu)準(zhǔn)確度的平均值。將Pi和Ci的加權(quán)組合作為啟發(fā)式信息，修改w的值可以調(diào)整Pi和Ci權(quán)重。

3.1.6 蟻群算法構(gòu)造信息SNP子集的偽代碼

以下使用偽代碼來描述用于選SNP 子集的蟻群算法的總體框架。當(dāng)方法達(dá)到最大迭代次數(shù)或達(dá)到重構(gòu)準(zhǔn)確度時退出。

基于蟻群算法的信息SNP 子集構(gòu)造

輸入：個體基因型數(shù)據(jù)或單倍型數(shù)據(jù)

輸出：信息SNP 集合

Begin：

初始化信息素和參數(shù);

Nc=0；

While(Nc＜=Nc_max)

For i=1 to n_ants //n_ants為螞蟻數(shù)量

//每只螞蟻分別逐個添加SNP到候選信息SNP

While(prediction accuracy is not enough)

計算ηi；

按式（4）在候選SNP選擇新的SNP位點(diǎn);

End while

保存本次迭代過程中的最優(yōu)解；

End For

2.1 兩組血漿NT-proBNP水平比較 病例組患兒的NT-proBNP水平在治療3 d、治療7d 、治療14 d時均顯著低于組內(nèi)治療前(F=176.405，P<0.05)，病例組患兒的NT-proBNP水平在治療前、治療3 d、治療7 d、治療14 d時均顯著高于對照組(F=286.557，P<0.05)。見表1。

//每只螞蟻結(jié)束尋找路徑后，根據(jù)已經(jīng)經(jīng)過的節(jié)點(diǎn)數(shù)目

計算這些節(jié)點(diǎn)上每只螞蟻留下的信息素

按式（6）計算新的信息素;

Nc++；

End while

返回候選SNP;

End

3.2 基于KNN的SNP預(yù)測

本文采用K-最近鄰（KNN）方法來預(yù)測未選擇的SNP。給定測試集d（其類別未知），該方法是在訓(xùn)練集中查找k 個最近鄰居，并且使用k 個最近鄰居的類別來對候選者進(jìn)行預(yù)測。在本文中，兩個SNP 基因序列之間的距離為漢明距離，必須確定基因型之間k 個最近的鄰居。本文設(shè)定為5-NN（k=5），即確定投票過程的基因型樣本的5 個鄰居，并通過對這五個鄰居進(jìn)行投票獲得預(yù)測樣本結(jié)果。整體過程如圖1所示。

4 實驗

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

圖1 SNP選擇流程圖

實驗數(shù)據(jù)本實驗中所使用的數(shù)據(jù)由無錫市精神衛(wèi)生中心提供。數(shù)據(jù)格式為遺傳病的SNP 基因型格式，并且每個樣本都帶有標(biāo)記信息，標(biāo)記樣本是否患病。數(shù)據(jù)集的概況描述如表1所示。

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

4.2 實驗評價指標(biāo)

本文使用信息SNP 子集對非信息SNP 子集的重構(gòu)準(zhǔn)確度（ACC(I)）作為信息SNP子集的評價指標(biāo)，其定義如式（10）所示：

其中g(shù)i為位點(diǎn)上的實際值，為預(yù)測出的值，兩者差的絕對值即為預(yù)測誤差。N 是樣本的數(shù)量，是非信息SNP 的數(shù)量，兩者積為總重構(gòu)次數(shù)。重構(gòu)度越高，信息SNP子集對非信息SNP的預(yù)測效果越好。

4.3 SNP數(shù)據(jù)編碼

由于實驗的原始數(shù)據(jù)是SNP的基因型表示，本次實驗采用的編碼方式是“0-1-2”編碼，分別表示AA、Aa以及aa。編碼完成后，對缺失值需要進(jìn)一步進(jìn)行填充?？紤]到SNP局部可能存在關(guān)聯(lián)性，所以使用K近鄰的方式對其進(jìn)行填充。

4.4 實驗結(jié)果與分析

Halperin 等設(shè)計了一種最大投票法STAMPA，將此方法用于對結(jié)果的重構(gòu)，每個位點(diǎn)根據(jù)其最相似的信息SNP 位點(diǎn)預(yù)測［18］。粒子群算法BPSO［19］類似于本文中候選子集構(gòu)造的蟻群算法，BPSO 和STAMPA 組合為BPSO/STAMPA，與MLR 的組合是BPSO/MLR，本文方法蟻群算法與最近鄰分類組合標(biāo)記為ACO/KNN，蟻群算法參數(shù)設(shè)置為α=1、β=3，螞蟻數(shù)量設(shè)置為10，迭代次數(shù)為10。

4.4.1 重構(gòu)準(zhǔn)確度

三種方法在兩個數(shù)據(jù)集上的實驗結(jié)果如圖2和圖3 所示。在圖中，橫坐標(biāo)是信息SNP 的數(shù)量，縱坐標(biāo)是重構(gòu)準(zhǔn)確率?？梢钥闯觯谠摂?shù)據(jù)集E144 中，所提出的方法和BPSO/MLR 的準(zhǔn)確度明顯優(yōu)于BPSO/ STAMPA，ACO/KNN 與BPSO/MLR的重構(gòu)準(zhǔn)確度效果相似。在數(shù)據(jù)集G1000中，該方法具有比BPSO/STAMPA 和BPSO/MLR 的更高的重構(gòu)準(zhǔn)確度，并且重構(gòu)準(zhǔn)確度平均高出2%～5%。BP?SO/MLR 略高于BPSO/STAMPA，ACO/KNN 引入了連鎖度，使得位點(diǎn)間的連鎖不平衡性較高。可以得出本文方法更能獲取有利于樣本重構(gòu)的特征。

圖2 數(shù)據(jù)集E144上重構(gòu)準(zhǔn)確度

4.4.2 運(yùn)行時間

分別使用ACO/KNN、BPSO/STAMPA 和BPSO/MLR 三種算法在兩個數(shù)據(jù)集上進(jìn)行實驗，并比較每種算法的運(yùn)行時間，其結(jié)果如圖4和圖5所示（每組實驗重復(fù)三次并對結(jié)果取均值）。在圖中，橫坐標(biāo)是信息SNP 的數(shù)量，縱坐標(biāo)是運(yùn)行時間（單位：s）。由于E144 數(shù)據(jù)集規(guī)模大于G1000 數(shù)據(jù)集，因此運(yùn)行時間也大于G1000 數(shù)據(jù)集。從圖中可以看出，隨著SNP 數(shù)量的增加，該方法的優(yōu)勢也逐漸明顯。

圖3 數(shù)據(jù)集G1000上重構(gòu)準(zhǔn)確度

圖4 數(shù)據(jù)集E144上運(yùn)行時間

圖5 數(shù)據(jù)集G1000上運(yùn)行時間

5 結(jié)語

在本文中，針對SNP 數(shù)據(jù)普遍存在的少樣本、高維度的問題，和不同SNP位點(diǎn)之間存在連鎖不平衡導(dǎo)致的位點(diǎn)之間具有強(qiáng)相關(guān)性的特點(diǎn)，將連鎖不平衡性應(yīng)用到蟻群算法中，提出一種基于蟻群算法（ACO）信息SNP 選擇方法。本文使用的實驗數(shù)據(jù)來自無錫市精神衛(wèi)生中心，并與文獻(xiàn)中的SNP選擇方法作了比較。本文的后續(xù)工作是對KNN 進(jìn)行改進(jìn)，使篩選出的信息SNP子集具有更高的重構(gòu)度。