吳蓉暉，盧友敏，2?

（1.湖南大學(xué) 信息科學(xué)與工程學(xué)院，湖南 長沙 410082；2.懷化學(xué)院 計算機工程系，湖南 懷化 418400）

隨著人類基因組計劃（HGP）測序工作的完成，生命科學(xué)的研究重點已經(jīng)從確定DNA序列組成轉(zhuǎn)移到了研究基因功能.由于復(fù)雜疾病［1］在人群中具有高死亡率及難以治愈等特點，使得復(fù)雜疾病成為醫(yī)學(xué)、生物學(xué)相關(guān)科研人員的重點研究對象.復(fù)雜疾病不同于孟德爾疾病，它的形成與發(fā)展通常涉及到多個基因的相互作用或者基因與環(huán)境的交互作用即上位作用.而從分子層次上看，上位作用即為基因調(diào)控網(wǎng)絡(luò)或生物化學(xué)代謝通路中的生物分子（例如DNA，RNA的蛋白質(zhì)等）之間的物理相互作用［2］.通常上位作用在基因型和疾病表型之間一般都表現(xiàn)為非線性關(guān)系，從而難以被檢測.在特殊情況下，單個基因與表型之間并沒有表現(xiàn)出相關(guān)性，但是當(dāng)該基因與其他基因或者環(huán)境聯(lián)合分析時，則存在明顯的上位作用.因此，復(fù)雜疾病一般具有表型異質(zhì)性、遺傳異質(zhì)性等特點，使人們難以從根本上理解其致病機理.

全基因組范圍內(nèi)的復(fù)雜疾病易感基因的發(fā)掘及其與疾病關(guān)聯(lián)方式的確定，將有利于更全面地理解復(fù)雜疾病發(fā)病機理，從而實現(xiàn)復(fù)雜疾病的預(yù)防、診斷和治療.盡管針對復(fù)雜疾病的SNP芯片已經(jīng)產(chǎn)生海量的數(shù)據(jù)，但是由于該數(shù)據(jù)本身具有的特征維數(shù)高和上位性分析存在組合爆炸等特點，使得該研究中如何對數(shù)據(jù)進行有效降維，并保留關(guān)鍵的上位作用，并有效刻畫上位作用與復(fù)雜疾病之間關(guān)系，成為了復(fù)雜疾病的全基因組關(guān)聯(lián)研究的熱點.本文首先采用多準(zhǔn)則融合策略對無關(guān)、冗余SNP位點進行過濾，然后采用蟻群優(yōu)化算法進一步剔除冗余SNP位點，實現(xiàn)對數(shù)據(jù)的降維并找出與疾病相關(guān)的上位性組合，然后采用支持向量機作為分類模型.實驗表明，本文方法具有實用意義.

1 全基因組關(guān)聯(lián)研究中存在的問題

復(fù)雜疾病上位性研究一般由數(shù)據(jù)預(yù)處理，上位性檢測以及分類評估3個階段構(gòu)成.目前，研究人員在這3個階段，提出將關(guān)聯(lián)統(tǒng)計分析、機器學(xué)習(xí)等方法應(yīng)用到該研究中，從而發(fā)展了很多上位性分析的模型及算法.

圖1 基于蟻群算法的分析框架Fig.1 The analysis framework

統(tǒng)計檢驗方法［3］如信息增益、方差檢驗和卡方檢驗等被用于上位性檢測，這些方法都暗含了各個特征SNP之間是相互獨立的這一假設(shè)，因此，在對特征進行評估時只考慮了特征與疾病性狀之間的關(guān)系，而忽略了特征與特征之間的相互作用，因此對于評估可能包含上位作用的生物數(shù)據(jù)存在準(zhǔn)確率較低等缺點.在眾多機器學(xué)習(xí)算法中，決策樹算法是較早被確認為是識別SNP－SNP相互作用的有用工具，但是該類方法只應(yīng)用于相對較小的數(shù)據(jù)集.為了解決較大規(guī)模數(shù)據(jù)集上的上位性分析，Chen等［4］研究了隨機森林中的統(tǒng)計效率，用于分析包含了成百上千個候選SNPs的疾病數(shù)據(jù)集.目前，雖然這些方法具有一定的優(yōu)勢，但是仍然存在時間復(fù)雜度高、分類準(zhǔn)確率低、假陽性高等不足.針對以上存在的問題，當(dāng)前已有一些研究提出先過濾掉冗余、無關(guān)的SNP再進行上位性分析的策略.如果單個SNP對疾病具有統(tǒng)計可檢測的主效應(yīng)，那么可以檢測出其與疾病之間的關(guān)聯(lián)（association），然后過濾掉低關(guān)聯(lián)強度的SNP，從而縮小后續(xù)上位性組合檢測中所需搜索的組合空間.但是，某種情況下可能出現(xiàn)純上位性現(xiàn)象，傳統(tǒng)單SNP分析方法可能剔除了這些位點，而導(dǎo)致后續(xù)分類準(zhǔn)確低，并且由于組合爆炸，對所有SNP組合進行窮舉搜索大大增加了計算復(fù)雜度.因此，亟需一種有效的篩選方法剔除無關(guān)、冗余SNP位點，并有效保留純上位作用SNP位點及主效應(yīng)SNP位點，在保證分類準(zhǔn)確度的基礎(chǔ)上降低計算復(fù)雜度.

2 基于蟻群算法的分析方法

其中S′表示構(gòu)成上位性的SNP組合；C（S′）為SNP組合G′的分類能力.在具有相同分類能力的不同SNP組合之間選擇較小的SNP組合，是符合復(fù)雜疾病研究的發(fā)展規(guī)律.

由于復(fù)雜疾病SNP芯片數(shù)據(jù)具有的特征維數(shù)高以及上位作用等特點，為了保證在對SNP數(shù)據(jù)進行降維基礎(chǔ)上，同時保留其中的關(guān)鍵上位組合，并有效對SNP上位性組合與復(fù)雜疾病之間建立映射，本文提出先過濾后分類的分析框架，如圖1所示.

基于分類的復(fù)雜疾病上位性分析有一個基本假設(shè)：如果某一SNP上位性組合與復(fù)雜疾病的形成相關(guān)，那么通過分析個體在這些SNP位點上的組合模式，則可以判別個體是否患病.利用該假設(shè)，對復(fù)雜疾病上位性分析轉(zhuǎn)換為以下數(shù)學(xué)模型：

2.1 多準(zhǔn)則融合過濾

在SNP數(shù)據(jù)的上位性分析中，面臨的最大挑戰(zhàn)是SNP組合空間的爆炸，而對所有的SNP組合進行窮舉分析，則是NP難問題.為了降低上位性分析時間復(fù)雜度，一種有效策略是對復(fù)雜疾病的患?。瓕φ諗?shù)據(jù)進行分析，然后利用某種過濾規(guī)則去除噪聲、無關(guān)SNP位點，這些位點主要表現(xiàn)為在對照樣本中與患病樣本中的SNP基因型基本一致，則可以認為它們是與復(fù)雜疾病無關(guān)的位點.

為了防止一些易感SNP因為單位點的弱效應(yīng)被過濾，本文提出采用多準(zhǔn)則融合策略綜合地、更為全面地評價每個SNP位點.主要原因有兩點：第一，借鑒集成多個弱分類器可以顯著提高分類的能力這一事實，采用多準(zhǔn)則融合可以更為準(zhǔn)確地評價每個SNP組合，從而降低假陽性；第二，因每個評價規(guī)則都具有獨特的傾向性，從而導(dǎo)致容易陷入局部最優(yōu)，而通過融合多種特征，可以更好地尋找全局最優(yōu)的上位性組合.本文采用對信噪比［5］、Relief［6］和卡方檢驗［7］準(zhǔn)則進行融合的方法來對SNP數(shù)據(jù)篩選過濾.

1）信噪比.廣義來講，信噪比（Signal to Noise Ratio）是指有效信息被破壞的程度，本文中用該指標(biāo)作為度量每個SNP位點對樣本分類貢獻的大小.

式中：d（s）為SNP位點s的打分值；μ＋s和μ－s分別為不同類別中s的基因型平均值；σ＋s和σ－s為基因型的標(biāo)準(zhǔn)差.從式（2）可以看出，SNP位點打分值越高，表明它在不同類間的差異越大或類內(nèi)變化率越小，那么其對于分類的貢獻越大.

2）Relief可作為一種基于權(quán)值的單位點排序方法，它通過多次迭代來評價位點的相關(guān)性，每次迭代過程中，首先隨機地從數(shù)據(jù)集中選擇一個樣本X，以及同一類中的X與它最近的鄰居H 和不同類中與X最近的M，然后利用公式（3）計算H 與M 的差別，從而更新所有特征的相關(guān)性.

式中：W（j）為SNP位點j在X 與目標(biāo)之間的相關(guān)性，迭代初始時被設(shè)置為0；dif f（j，x，x′）為SNP位點j在樣本x與x′上的差別.

3）Pearson卡方檢驗是Karl Pearson提出的用于檢驗樣本中某一些事件發(fā)生的概率是否等于理論分布的一種檢驗方法，也稱為擬合優(yōu)度檢驗.對于一般的I×J列聯(lián)表，Pearson檢驗可以表示為：

式中：Oi是第i個格子基因型觀測到的頻數(shù)；Ei是該基因型的理論頻數(shù)；∑對所有的格子求和.

由于以上3個標(biāo)準(zhǔn)對每個SNP位點評價的度量值不統(tǒng)一，本文首先對所有SNP位點分別按照以上3種不同準(zhǔn)則所對應(yīng)的重要程度進行排序，則每個位點i將有3種排名值分別為di，wi以及xi，然后將排名值相加得一個Si并排序，則并排序的名次綜合反映了每個SNP位點的重要程度，然后設(shè)定一個閾值（本文設(shè)為100），大于該閾值的位點則被過濾掉.

2.2 上位性分析

蟻群算法［8］已成功應(yīng)用于各個領(lǐng)域中的組合優(yōu)化NP難問題，如旅行商問題、圖著色問題以及微陣列特征選擇等.它具有天然的并行性，通過并行策略能極大提高運算速度.

人工蟻群算法由多個并行的螞蟻構(gòu)成，螞蟻之間通過概率密度函數(shù)進行通信，該函數(shù)由權(quán)重因子以及信息素濃度構(gòu)成.在蟻群算法用于上位性分析中，第k次迭代中位點i的選擇概率被定義為：

式中：α為信息素權(quán)值；β為啟發(fā)因子的權(quán)值；τi為第i個位點上信息素濃度；每個位點的啟發(fā)性信息ηi都被置為常數(shù)1.初始化時，每個位點都設(shè)置為相等的初始濃度值τ0.利用式（6），每只螞蟻m從所有SNP位點中選擇n（1≤n≤SNP位點數(shù)－1）只螞蟻分別構(gòu)造一個SNP上位性組合Sm，而每個SNP組合的分類性能則作為下一輪迭代中信息素更新的依據(jù)，更新函數(shù)為式（7）.其中，該過程中采用支持向量機［9］作為分類學(xué)習(xí)模型.

式中：ρ為大于0小于1的信息素揮發(fā)因子；Δτi（t－1）為第t－1次迭代中的最佳上位性組合Smax的分類準(zhǔn)確率.如果位點i屬于Smax，那么在第t次迭代中則按照式（7）改變其信息素，如果不屬于，則Δτi（t－1）等于0.以下是蟻群算法的偽代碼.

輸入：復(fù)雜疾病SNP數(shù)據(jù)集

輸出：SNP上位性組合位點

Step1：數(shù)據(jù)預(yù)處理；

Step2：初始化蟻群算法參數(shù)如蟻群規(guī)模iAntCount，最大迭代次數(shù)maxIteration；

Step3：每只螞蟻m根據(jù)概率選擇函數(shù)構(gòu)造上位性組合Sm；

Step4：利用支持向量機對每個Sm采用五折交叉驗證法評價分類性能；

Step5：記錄最優(yōu)分類性能的上位性組合Smax，更新每個位點信息素；

Step6：判斷是否滿足終止條件，如果不滿足則回到Step3，否則執(zhí)行Step7；

Step7：輸出上位性組合位點，退出程序.

2.3 分類評估方法

支 持 向 量 機 （Support Vector Machines，SVM）［9］是一種成熟的模式識別模型，它遵循結(jié)構(gòu)風(fēng)險最小化原則，在小樣本學(xué)習(xí)中體現(xiàn)出卓越優(yōu)勢，并且，其計算復(fù)雜度僅僅與支持向量數(shù)目有關(guān)，而與輸入空間維數(shù)無關(guān)，因此，它非常適宜處理復(fù)雜疾病SNP芯片數(shù)據(jù)這種典型的高維、少樣本數(shù)據(jù).

為了降低峰值現(xiàn)象，從而更可信地度量分類準(zhǔn)確率，本文采用五折交叉驗證法.五折交叉驗證法首先將樣本數(shù)據(jù)集分為5個子集，然后將其中4個作為訓(xùn)練集，另外一個作為測試集，進行一次分類測試，每個子集將被用作一次測試集，依次循環(huán)迭代5次，最后對5次分類準(zhǔn)確度求平均值，以此評價上位性SNP組合的分類準(zhǔn)確度.平均值的計算方法為：

式中：pi為第i次迭代的分類準(zhǔn)確度.

3 仿真實驗及分析

為了合理地評價該改進方法在上位性分析中的有效性，分別在分類準(zhǔn)確率以及運行時間指標(biāo)上對本文方法進行了驗證評價.首先采用C＋＋實現(xiàn)了本文算法，然后在 WIN7環(huán)境下執(zhí)行測試，測試環(huán)境的硬件配置為2G內(nèi)存，AMD雙核2.80GHz.

3.1 數(shù)據(jù)集

由于有采用的成本及涉及患者隱私的情況存在，使得復(fù)雜疾病分析的真實數(shù)據(jù)集中樣本量小，同時，有些真實復(fù)雜疾病數(shù)據(jù)中真正的致病基因，或者不同實驗分析的結(jié)果存在不一致性，因此無法用于驗證生物信息學(xué)方法得到的結(jié)果，所以現(xiàn)有研究中通常采用仿真數(shù)據(jù)評價機器學(xué)習(xí)方法.仿真數(shù)據(jù)生成需要設(shè)置幾個重要參數(shù)，如基因的外顯率函數(shù)、遺傳度以及次要等位基因頻率（MAF）等.本文仿真數(shù)據(jù)集來自參考文獻［10］，下載地址為 http：／／discovery.dartmouth.edu／epistatic＿data／.本文下載了2種不同參數(shù)的數(shù)據(jù)集，數(shù)據(jù)集的詳細介紹列于表1.

表1 模擬數(shù)據(jù)集Tab.1 Simulated datasets

3.2 實驗分析

利用以上數(shù)據(jù)集，用本文方法與SNPRuler算法［11］進行比較驗證.SNPRuler算法是基于預(yù)測規(guī)則推理和兩階段（two－stage）策略設(shè)計的，通過預(yù)測規(guī)則學(xué)習(xí)特征與類變量之間的關(guān)系，然后在測試數(shù)據(jù)上預(yù)測類標(biāo)簽.其上位性檢測中利用規(guī)則學(xué)習(xí)，原因在于：首先，上位性組合蘊含了一些模式或預(yù)測規(guī)則；再者，評估規(guī)則的尋找更為容易，更快捷.因此，SNPRuler方法通過挖掘預(yù)測規(guī)則來發(fā)現(xiàn)潛在的上位性組合.

SNP組合的上位性可以通過個體性狀的分類性能來評價，本文采用五折交叉驗證法驗證不同SNP組合的分類準(zhǔn)確度.在以上2個數(shù)據(jù)集中的分類準(zhǔn)確度的實驗結(jié)果分別如圖2和圖3所示.從圖2和圖3可知，在不同數(shù)據(jù)集上，本文方法的分類準(zhǔn)確度平均高于SNPRuler算法2%.通過分析可以發(fā)現(xiàn)，外顯率對上位性分析也有影響，外顯率高則更容易發(fā)現(xiàn)真正的致病位點，分類準(zhǔn)確度更高.

圖2 Dataset1的分類準(zhǔn)確度比較Fig.2 The prediction accuracy on Dataset1

圖3 Dataset2的分類準(zhǔn)確度比較Fig.3 The prediction accuracy on Dataset2

2種方法運行時間比較結(jié)果如圖4所示，運行時間分別對應(yīng)著不同方法尋找到最優(yōu)上位性組合即具有最高分類準(zhǔn)確度SNP組合所消耗的時間.由圖4可知，本文方法通過先過濾掉大量無關(guān)SNP后再搜索，使得上位性分析的運行時間總體看來大致接近SNPRuler算法的一半，較大地提高了上位性分析的效率.

圖4 運行時間比較Fig.4 The comparison of running time

4 結(jié)束語

為了探索與復(fù)雜疾病發(fā)生、發(fā)展相關(guān)的上位性，針對現(xiàn)有上位性分析方法存在高運算成本、假陽性高等不足，本研究提出了一種基于蟻群算法的上位性分析方法，它包含過濾以及上位性分析兩個階段，在過濾階段剔除大量無關(guān)位點后，使上位性分析過程的SNP組合空間大大縮小，使得高階上位性分析成為可能.并且，過濾階段采用了多準(zhǔn)則融合策略，更為全面、綜合地評價每個SNP位點，能有效保留單個弱效SNP位點.實驗表明，本文方法在分類準(zhǔn)確度以及運行時間上都有一定程度提高，具有實用意義.

［1］孫玉琳，趙曉航.復(fù)雜疾病基因定位策略與腫瘤易感基因鑒定［J］.生物化學(xué)與生物物理進展，2005，32（9）：804－809.SUN Yu－lin，ZHAO Xiao－h(huán)ang.The genetic mapping of complex diseases and the identification of tumors’susceptible genes［J］.Prog Biochem Biophys，2005，32（9）：804－809.（In Chinese）

［2］王文菊，尹先勇，崔勇，等.IL23／Th17通路基因上位性作用與漢族人銀屑病易感性研究［J］.實用醫(yī)院臨床雜志，2013，10（1）：1－3.WANG Wen－ju，YIN Xian－yong，CUI Yong，et al.Genes in IL23／Th17pathway have epistatic effects on psoriasis susceptibility in Chinese Han population［J］.Practical Journal of Clinical Medicine，2013，10（1）：1－3.（In Chinese）

［3］GENIN E，COUSTET B，ALLANORE Y，et al.Epistatic Interaction between BANK1and BLK in rheumatoid arthritis：results from a large trans－ethnic meta－analysis［J］.Plos One，2013，8（4）：e61044.

［4］CHEN S H，SUN J，DIMITROV L，et al.A support vector machine approach for detecting gene－gene interaction［J］.Genetic Epidemiology，2008，32：152－167.

［5］阮曉鋼，晁浩.腫瘤識別過程中特征基因的選?。跩］.控制工程，2007，14（4）：374－375.RUAN Xiao－gang，CHAO Hao.Selection of feature genes in cancer classification［J］.Control Engineering of China，2007，14（4）：374－375.（In Chinese）

［6］ROBNIK－?IKONJA M，KONONENKO I.Theoretical and empirical analysis of relief and relief［J］.Machine Learning，2003，53（1）：23－69.

［7］WAN Xiang，YANG Can，YANG Qiang，et al.The complete compositional epistasis detection in genome－wide association studies［J］.BMC Genetics，2013，14（7）：1－11.

［8］吳建輝，章兢，劉朝華.基于蟻群算法和免疫算法融合的TSP問題求解［J］.湖南大學(xué)學(xué)報：自然科學(xué)版，2009，36（10）：82－85.WU Jian－h(huán)ui，ZHANG Jin，LIU Zhao－h(huán)ua.Solution of TSP problem based on the combination of ant colony algorithm and immune algorithm［J］.Journal of Hunan University：Natural Sciences，2009，36（10）：82－85.（In Chinese）

［9］文益民，王耀南，張瑩.基于分類面拼接的快速模塊化支持向量機研究［J］.湖南大學(xué)學(xué)報：自然科學(xué)版，2009，36（3）：46－49.WEN Yi－min，WANG Yao－nan，ZHANG Ying.On pasting small fast modular SVMs for classification［J］.Journal of Hunan University：Natural Sciences，2009，36（3）：46－49.（In Chinese）

［10］WANG Y，LIU G M.An empirical comparison of several recent epistatic interaction detection methods［J］.Bioinformatics，2011，27（21）：2936－2943.

［11］WAN Xiang，YANG Can，YANG Qiang，et al.Predictive rule inference for epistatic interaction detection in genome－wide association studies［J］.Bioinformatics，2010，26（1）：30－37.