一種面向大規(guī)模序列數(shù)據(jù)的交互特征并行挖掘算法

趙宇海 印 瑩 李 源 汪嗣堯 王國仁

(東北大學(xué)計算機(jī)科學(xué)與工程學(xué)院 沈陽 110819)

序列數(shù)據(jù)是一種重要的數(shù)據(jù)類型，被廣泛應(yīng)用于科學(xué)與工程學(xué)、商業(yè)、客戶行為分析、股票趨勢預(yù)測、DNA序列分析、Web使用行為分析等諸多實際應(yīng)用領(lǐng)域[1-2].大數(shù)據(jù)時代，信息技術(shù)的不斷廉價化與互聯(lián)網(wǎng)及其延伸所帶來的無處不在的信息技術(shù)應(yīng)用，進(jìn)一步加劇了海量序列數(shù)據(jù)的積累.以生物領(lǐng)域為例，高通量測序技術(shù)，尤其是新一代測序技術(shù)的迅猛發(fā)展，生物序列數(shù)據(jù)的規(guī)模呈爆炸性增長.大規(guī)模生物序列分析成為當(dāng)前生物信息領(lǐng)域的熱點研究問題之一[3-5].

基于單核苷酸多態(tài)性(single nucleotide polymor-phism, SNP)數(shù)據(jù)的復(fù)雜疾病全基因組關(guān)聯(lián)研究是大規(guī)模生物序列分析的典型應(yīng)用案例之一，其目的是在基因組范圍內(nèi)以SNP為標(biāo)記，通過患病人群和正常人群之間的對比，找出顯著序列變異，然后從中篩選出與疾病相關(guān)的位點(locus).該問題本質(zhì)上可建模為特征選擇問題，即從數(shù)據(jù)中選擇出與分類屬性關(guān)聯(lián)性最強(qiáng)的特征子集.雖然問題簡單明了，但解決起來卻并不輕松，面臨著一些嚴(yán)峻的挑戰(zhàn).當(dāng)前發(fā)現(xiàn)的致病SNP中只有少數(shù)能中等或大量地增加復(fù)雜疾病的致病風(fēng)險，一些已經(jīng)被生物實驗證實的致病SNP位點locus沒有被GWAS(genome-wide associa-tion study)識別出來.研究結(jié)果僅解釋了疾病很小部分的遺傳變異，出現(xiàn)了“遺傳丟失”現(xiàn)象，甚至彼此矛盾[3,6].例如，目前已發(fā)現(xiàn)了超過40個位點的變異影響人類身高性狀，但在幾萬個樣本實驗中發(fā)現(xiàn)這些位點大概僅解釋了5%左右的表型變異[3]；類似地，在糖尿病研究中發(fā)現(xiàn)的18個位點僅解釋了6%左右的遺傳風(fēng)險；再以兒茶酚氧位甲基轉(zhuǎn)移酶基因COMT(catechol-O-methyltransferase)為例，Kotler等人[7]的研究支持158A(Met)會升高精神分裂癥的發(fā)病率，而Wonodi等人[8]的研究則支持它會降低發(fā)病率.

造成以上情況的主要原因之一在于，傳統(tǒng)GWAS研究僅僅關(guān)注了單點關(guān)聯(lián)分析，忽略了微效位點間的交互作用對疾病的影響.生物體是個復(fù)雜系統(tǒng)，內(nèi)部元件之間存在著錯綜復(fù)雜的關(guān)系，研究位點間的交互作用是解析復(fù)雜疾病性狀遺傳機(jī)理的重要途徑之一.很多遺傳學(xué)家相信位點間的交互作用在復(fù)雜疾病性狀的遺傳變異中扮演著重要角色.越來越多的理論和實驗證據(jù)表明：SNP交互作用是復(fù)雜疾病產(chǎn)生和發(fā)展的重要遺傳學(xué)基礎(chǔ)之一[3-5].研究者們已在對老年癡呆癥、乳腺癌、精神分裂癥等復(fù)雜疾病的研究中發(fā)現(xiàn)了真實的SNP交互作用.全基因組SNP交互作用研究為GWAS開辟了嶄新的研究思路和科研方向[1,5,9].針對復(fù)雜疾病，開展基于交互作用的高通量生物數(shù)據(jù)挖掘與分析技術(shù)研究成為新一波生物信息研究浪潮中人們關(guān)注的焦點[1,5,9].

SNP交互作用發(fā)現(xiàn)面臨著嚴(yán)峻的計算挑戰(zhàn).全基因組關(guān)聯(lián)研究中的SNP數(shù)目動輒幾十萬到上百萬.即使只檢驗2位點互作，在30萬個位點上的檢測次數(shù)也已高達(dá)1 000億次，對3階甚至更高階的交互作用更是難以想象的天文數(shù)字.例如假設(shè)每秒可以檢測105個SNP位點組合，則僅3階交互作用計算就需要約1.4×108年.因此，雖然目前生物上已證實高階交互作用在復(fù)雜疾病中真實存在[10-11]，但鮮有涉及該領(lǐng)域的研究報道.

針對大規(guī)模生物序列中高階交互特征挖掘面臨的密集計算問題，本文提出了一種新的基于并行處理和演化計算的解決框架，大致由4個主要步驟組成：

1) 數(shù)據(jù)預(yù)處理. 提出一種基于塊和位點的2階段過濾方法，提前過濾無關(guān)特征，提高整個框架的執(zhí)行效率.

2) “低耦合高內(nèi)聚”的數(shù)據(jù)劃分. 執(zhí)行序列挖掘任務(wù)，將原始高維特征數(shù)據(jù)劃分成一系列低維特征區(qū)域，并保證其相互間的低交互可能性(即每個區(qū)域內(nèi)的特征交互可能性高，不同區(qū)域間的特征交互可能性低).

3) 多樣性的特征區(qū)域篩選. 許多劃分后的特征區(qū)域具有較高相似度，如果在所有特征區(qū)域上都執(zhí)行交互特征挖掘，會產(chǎn)生很多重復(fù)結(jié)果，降低計算效率. 將采用最小獨立支配集的思想，基于多樣性原則進(jìn)行代表性特征區(qū)域選擇.

4) 交互特征挖掘.提出基于置換搜索的并行蟻群算法，進(jìn)行高效且有效的交互特征選擇.

本文的主要貢獻(xiàn)有4個方面：

1) 提出了一種面向大規(guī)模序列數(shù)據(jù)的交互特征并行挖掘框架.現(xiàn)有框架大多從減少各節(jié)點處理的序列數(shù)角度提高計算效率，因此對數(shù)據(jù)進(jìn)行水平分組.與此不同，提出的框架從減少分組數(shù)據(jù)位點數(shù)入手，對數(shù)據(jù)進(jìn)行垂直分解.位點數(shù)是制約互作挖掘效率的根本因素，提出的框架具有天然的優(yōu)勢.

2) 提出了極大等位公共子序列(maximal allelic common subsequence, MACS)的概念并設(shè)計了基于MACS的特征區(qū)域劃分策略.該策略能保證交互特征在劃分區(qū)域間的低耦合性和劃分區(qū)域內(nèi)的高內(nèi)聚性，避免計算中的耦合引起的高額通信代價.

3) 設(shè)計并實現(xiàn)了一種基于置換策略的并行蟻群局部搜索算法(parallel ant colony optimization for interactive features search, PACOIFS)，避免了特征子集搜索過程中的邊際作用影響并保持計算的高效性.

4) 通過將本文算法和其他3個算法在糖尿病患者數(shù)據(jù)集和模擬數(shù)據(jù)集上進(jìn)行測試分析，證明了算法的高效性和生物有效性.

1 相關(guān)工作

根據(jù)搜索策略，目前的交互特征挖掘方法可分為：枚舉搜索、隨機(jī)搜索和貪婪搜索方法.枚舉搜素對所有特征組合進(jìn)行檢測，優(yōu)點是結(jié)果準(zhǔn)確，但計算量大，不適用于大規(guī)模數(shù)據(jù)；隨機(jī)搜索隨機(jī)地從所有位點中挑選位點子集，考察當(dāng)前位點的重要性.重復(fù)以上步驟，盡可能覆蓋所有位點.當(dāng)位點總數(shù)大時，很難覆蓋所有位點；貪婪搜索算法不挑選盡可能覆蓋所有位點的多個位點子集，先通過某種標(biāo)準(zhǔn)過濾掉絕大部分噪音位點，再對剩余位點探索交互作用.

SNPHarvester算法[13]是一種基于啟發(fā)式搜索的算法，利用路徑選擇過程來采樣搜索空間.其基本思想是先通過單位點顯著性檢測刪除有主效應(yīng)的位點，在剩余位點中用爬山算法貪婪地搜索出一組互作效應(yīng)最強(qiáng)的位點集合，最后利用懲罰回歸從候選結(jié)果中提取交互作用.SNPHarvester計算速度快，每發(fā)現(xiàn)一個局部最優(yōu)后，對應(yīng)的位點集合就被從搜索空間中移除，搜索空間越來越小.其主要問題在于算法易陷入局部最優(yōu)且不易并行化.

AntEpiSeeker算法[14]是一種基于蟻群優(yōu)化的2階段啟發(fā)式搜索方法，可適用于大規(guī)模數(shù)據(jù)集中位點互作的檢測.在篩選階段，AntEpiSeeker用蟻群算法搜索含有多于給定位點數(shù)的SNPs子集，從中選出卡方統(tǒng)計量最大的若干子集作為候選子集，再選出信息素最大的若干位點作為候選位點.在候選子集與候選位點中，通過窮舉搜索發(fā)現(xiàn)包含給定數(shù)目位點的集合，采用卡方統(tǒng)計量的p-值作為評價函數(shù).以上2階段搜索重復(fù)2輪.AntEpiSeeker對邊際效應(yīng)魯棒性較好，但需要指定一系列參數(shù)(包括蟻群算法的參數(shù)和2輪迭代涉及的參數(shù))，算法性能受參數(shù)影響大.

2 基本概念與問題描述

2.1 基本概念

本節(jié)將給出文中所涉及到的一些基本概念并對要解決的問題給出形式化描述.

Fig. 1 An example of maximal allelic common subsequence圖1 極大等位公共子序列的示例圖

極大等位公共子序列與廣為熟知的最長公共子序列不同.如圖1所示，若S1=100101和S2=101011為任意2個由0和1構(gòu)成的序列，則序列Sα=101為S1和S2的極大等位公共子序列，序列Sβ=10101為S1和S2的最長公共子序列.顯然，Sα≠Sβ.計算多序列的最長公共子序列是典型的NP-難問題[1]，時間復(fù)雜度相對于序列數(shù)量N是指數(shù)級別的，而由圖1不難得知，極大等位公共子序列通過簡單的位“與”計算即可獲得，時間復(fù)雜度僅為O(N).因此，極大等位公共子序列的計算量遠(yuǎn)小于最長公共子序列.

由定義1可知，對給定的數(shù)據(jù)集D，其中所有任意k條序列的極大等位公共子序列的數(shù)量隨著序列數(shù)|D|的增加而增加.一些極大等位公共子序列可能會存在較大的序列相似性.如果考慮所有的極大等位公共子序列，將引入很多重復(fù)計算，影響交互特征挖掘的執(zhí)行效率.因此，提出基于序列多樣性的特征區(qū)域選擇方法.通過將問題規(guī)約到經(jīng)典的最小獨立支配集問題，設(shè)計有效的解決方法，使得在盡可能不丟失結(jié)果的前提下，減少很多重復(fù)計算.以下是最小獨立支配集的具體定義.

定義2.最小獨立支配集. 給定圖G={V,E}，其中，V和E分別表示該圖的點集和邊集.若存在V*?V，滿足?u,v∈V*，(u,x)?E,并且?w∈V,要么w∈V*，要么?x∈V*，(x,w)∈E，則稱V*為V的獨立支配集(independent dominating set, IDS).進(jìn)一步，包含節(jié)點數(shù)最少的獨立支配集被稱為最小獨立支配集(minimum independent dominating set, MIDS).

本工作的主要任務(wù)是從大規(guī)模序列中挖掘與目標(biāo)變量(如某種疾病或表型)存在顯著關(guān)聯(lián)的高階交互特征，因此有必要先明確k階交互特征的定義.

定義3.k階交互特征. 令F為含有k個特征f1,f2,…,fk的特征子集，C為某個特征或特征子集與目標(biāo)變量(類標(biāo)簽)相關(guān)性的度量.若滿足對F的任一劃分F={F1,F2,…,Fl}，C(F)>C(Fi)，其中i∈[1,l]，l≥2且Fi≠?，則稱f1,f2,…,fk為k階交互特征.

由定義3可知，如果f1,f2,…,fk為k階交互特征，當(dāng)且僅當(dāng)其對目標(biāo)變量(類標(biāo)簽)的影響大于其任一子集對目標(biāo)變量(類標(biāo)簽)的影響.

2.2 問題定義

統(tǒng)計顯著性是常用的特征質(zhì)量度量標(biāo)準(zhǔn).對給定的候選特征集Fi，基于單次假設(shè)檢驗即可判斷其與類標(biāo)簽cr是否存在顯著的統(tǒng)計關(guān)聯(lián).但是，大規(guī)模序列數(shù)據(jù)中的待檢驗k階交互特征數(shù)量極其龐大，多重假設(shè)檢驗引起的假陽性結(jié)果增多問題會給后續(xù)的生物驗證帶來巨大的資源浪費.族錯誤率FWER(family wise error rate)是指挖掘結(jié)果中至少存在一個假陽性結(jié)果的概率P(FP≥1)，其中FP代表總體結(jié)果中的假陽性個數(shù).FWER≤α意味著總體挖掘結(jié)果中至少出現(xiàn)一個假陽性結(jié)果的概率被控制在α以下.FWER在總體假陽性結(jié)果控制方面有堅實的理論基礎(chǔ)[15].因此，本研究要求挖掘結(jié)果滿足給定的族錯誤率閾值α.

問題描述：給定包含m條序列S={s1,s2,…,sm}和t個類標(biāo)簽C={c1,c2,…,ct}的數(shù)據(jù)集D，記si={fi 1,fi 2,…,fi n}為序列si包含的n個特征，cj為si對應(yīng)的類標(biāo)簽，其中1≤i≤m，1≤j≤t，發(fā)現(xiàn)：所有滿足定義3，與特定類標(biāo)簽cr∈C存在顯著統(tǒng)計關(guān)聯(lián)的k階交互特征并且滿足FWER≤α.

3 PACOIFS框架

第2.1節(jié)所述，大規(guī)模生物序列中的高階交互特征挖掘面臨著密集計算問題.本節(jié)將提出一種新的基于并行處理和演化計算的解決框架(PACOIFS)，并對構(gòu)成框架的各主要步驟加以詳細(xì)介紹.

3.1 整體思想

圖2給出了PACOIFS框架的總體示意圖:

Fig. 2 An illustration of PACOIFS圖2 PACOIFS框架示意圖

框架由4個主要步驟構(gòu)成：

1) 數(shù)據(jù)預(yù)處理.對原始數(shù)據(jù)進(jìn)行編碼后，執(zhí)行特征維度削減，提前過濾一些無關(guān)特征，提高整個框架的執(zhí)行效率.

2) 數(shù)據(jù)劃分.基于極大等位公共子序列MACS，將原始高維特征數(shù)據(jù)劃分成一系列低維特征區(qū)域，并保證其“低耦合高內(nèi)聚”性(即特征交互在區(qū)域內(nèi)的可能性高，在區(qū)域間的可能性低).

3) 特征區(qū)域篩選.許多劃分后的特征區(qū)域具有較高的相似度.如果在所有特征區(qū)域上都執(zhí)行交互特征挖掘，會產(chǎn)生很多重復(fù)結(jié)果，降低計算效率.采用最小獨立支配集的思想，基于多樣性原則選擇代表性的特征區(qū)域.

4) 交互特征挖掘.基于選定的代表模式對數(shù)據(jù)進(jìn)行垂直分解，提出基于置換搜索的并行蟻群算法，進(jìn)行高效且有效的交互特征選擇.

3.2 基于塊和位點過濾的數(shù)據(jù)預(yù)處理

原始SNP數(shù)據(jù)一般為2種形式：基因型或單體型數(shù)據(jù).以某位點為例，若用A表示出現(xiàn)頻率較高的等位基因，a表示出現(xiàn)頻率較低的等位基因，則基因型數(shù)據(jù)有3種形式：AA，Aa或aa，單體型數(shù)據(jù)有2種形式：A或a. 為便于計算機(jī)處理，通常需要先對原始SNP數(shù)據(jù)進(jìn)行編碼，前者的3種狀態(tài)分別編碼為0，1，2，后者的2種狀態(tài)分別編碼為0和1. 提出的方法可同時兼容基因型數(shù)據(jù)和單體型數(shù)據(jù).

理論上，對特定的復(fù)雜疾病，不可能所有的位點都與之具有遺傳學(xué)關(guān)聯(lián)[16]. 稱基因組中與特定疾病具有真正遺傳學(xué)關(guān)聯(lián)的區(qū)域為“熱區(qū)”.“熱區(qū)”定位是影響交互特征發(fā)現(xiàn)效率和有效性的關(guān)鍵因素.很多方法通過采用基于單位點過濾的方法，保留可能的“熱區(qū)”. 其問題在于，單位點過濾大多只保留具有強(qiáng)邊際效應(yīng)的候選位點，許多具有弱邊際效應(yīng)但具有強(qiáng)交互作用的位點組合被遺漏了.在對原始數(shù)據(jù)編碼后，本節(jié)中將提出一種結(jié)合塊過濾和位點過濾的2階段過濾方法BLFilter. 先通過基于圖論的塊過濾，粗粒度地保留與疾病具有較高關(guān)聯(lián)可能性的區(qū)域；再通過細(xì)粒度的位點過濾進(jìn)一步提煉保留區(qū)域內(nèi)的位點，以最大限度地保留邊際效應(yīng)弱但交互作用強(qiáng)的位點.

Fig. 3 The graph model of blocks圖3 分塊的圖模型

3.2.1 塊過濾

(1)

根據(jù)式(1)，圖G中的密集子圖對應(yīng)的區(qū)域分塊內(nèi)部和相互之間都存在較多的顯著交互位點對.因此，直觀上可理解為與疾病存在顯著關(guān)聯(lián)的“熱點”區(qū)域，問題被轉(zhuǎn)化為密集子圖發(fā)現(xiàn)問題.目前，從圖G中發(fā)現(xiàn)最密子圖的最優(yōu)算法時間復(fù)雜度為O(|V|(|V|+|E|)log(|V|)log(|V|+|E|))[17]，雖然是多項式級的，但對大規(guī)模序列分析而言，仍是不現(xiàn)實的.現(xiàn)有的2-近似算法貪婪算法[18]沒有考慮頂點和邊帶有權(quán)重的圖.針對該問題，本節(jié)提出一種能有效處理邊權(quán)重和點權(quán)重的密集子圖發(fā)現(xiàn)算法GREEDYVED.該算法迭代地移除當(dāng)前圖中具有最小平均度數(shù)的頂點vx及相關(guān)的邊，并計算移除后圖G′的密度.整個迭代結(jié)束后，d(G′)最大的子圖被作為最密子圖輸出.偽碼如算法1所示，其時間復(fù)雜度為O(|V|+|E|).

算法1.塊過濾算法GREEDYVED.

輸入：權(quán)重圖G=(V,E)；

①Sn=V;

② fori=nto 1 do

④Si-1=Si{vx}；

⑤ end for

定理1對塊過濾算法GREEDYVED的近似程度給出了證明.

(2)

證畢.

3.2.2 位點過濾

對保留區(qū)域內(nèi)的位點，進(jìn)一步執(zhí)行細(xì)粒度的位點過濾.偽碼如算法2所示:

算法2.位點過濾算法.

輸入：塊過濾后的位點集合X；

輸出：位點子集X′?X.

② fori=2 to |X| do

then

max_pos=i并且保留max_pos之前的所有位點至X′;

⑤ end if

⑥ end for

⑦ fori=max_pos+1 to |X| do

then

⑨ 將Xi和Xj保留至X′；

⑩ end if

具體過程如下：

與傳統(tǒng)方法相比，基于塊和位點的2階段過濾方法，既保持了過濾的高效性，又考慮了過濾的有效性.不僅保留了具有強(qiáng)邊際效應(yīng)的位點，而且最大限度的考慮了邊際效應(yīng)弱，但具有潛在交互作用的位點.

3.3 基于MACS的數(shù)據(jù)劃分

大規(guī)模序列中的高階交互特征挖掘是典型的計算密集型任務(wù)，位點數(shù)是制約交互作用挖掘效率的根本因素.現(xiàn)有并行框架大多對數(shù)據(jù)進(jìn)行水平分組，從減少各節(jié)點處理的序列數(shù)量角度提高計算效率.與此不同，本節(jié)從減少各節(jié)點處理的位點數(shù)入手，提出一種基于MACS的數(shù)據(jù)垂直分解策略，保證交互作用劃分后的“低耦合高內(nèi)聚”性，即每個區(qū)域內(nèi)的特征交互可能性高，不同區(qū)域間的特征交互可能性低.

3.3.1 劃分粒度

如問題定義所述，本研究要求挖掘結(jié)果的族錯誤率FWER被控制在α以下.利用統(tǒng)計中的置換檢驗方法，理論上可以獲得滿足FWER≤α的最大單次假設(shè)檢驗閾值δ，但計算量極其巨大.Llinares-López等人于2015年提出了一種基于假設(shè)檢驗的高效顯著模式挖掘算法FastWY[15]，一種快速確定δ的計算方法，并證明了模式的支持度γ與顯著性p值之間的對應(yīng)關(guān)系.

Table 1 p_value based on the Contingency表1 基于列聯(lián)表的p值計算

(3)

進(jìn)一步，證明了最小可達(dá)p值的下界ψ*(xi)相對于xi是單調(diào)遞減的.利用FastWY快速確定δ之后，若記其對應(yīng)的xi為ψ*-1(δ)，則只需檢測所有出現(xiàn)次數(shù)不低于ψ*-1(δ)的特征集合.根據(jù)定義1，數(shù)據(jù)集D中任意ψ*-1(δ)個序列的MACS構(gòu)成的子序列集合M即為劃分后的交互特征候選區(qū)間.注意：M中的任一序列s在D中的支持度至少為ψ*-1(δ)，支持度不低于ψ*-1(δ)的模式或者為M中的序列或者為M中序列的子序列.因此，顯著交互特征只可能存在于M中的每個序列內(nèi)部，不同序列之間的位點不可能存在顯著交互作用.

3.3.2 基于MACS的并行劃分

目前，任務(wù)規(guī)約為求數(shù)據(jù)集D中任意ψ*-1(δ)個序列的MACS構(gòu)成的子序列集合M的問題.本研究中利用基于MapReduce的并行框架來高效地求得結(jié)果.

3.4 基于MIDS的特征區(qū)域篩選

基于交集操作得到的MACS數(shù)量隨著數(shù)據(jù)行數(shù)的增多而顯著增多，在大規(guī)模序列數(shù)據(jù)分析時對整個算法的挖掘效率影響非常大.在實驗中發(fā)現(xiàn)了一個現(xiàn)象，即許多MACS具有較高相似性.如果不加處理，直接基于所有的MACS執(zhí)行挖掘過程，會產(chǎn)生大量冗余的計算.因此，在盡可能保證結(jié)果完整性的前提下對得到的特征區(qū)域進(jìn)行多樣性選擇.目的是選擇盡可能少的代表性區(qū)域，使其相互間差異大，但與未選擇的非代表性區(qū)域有較大的相似性.

該任務(wù)可以規(guī)約為最小獨立支配集(MIDS)問題.每個MACS作為圖G中的1個頂點，如果MACS間的相似度大于給定閾值，2個頂點間有邊存在.提出了基于貪心策略的(greedy,relevance,diversity,coverage, GRDC)算法進(jìn)行特征區(qū)域多樣性選擇.基本思想如下：每次從無向圖中選擇度最大的節(jié)點并將其加入最終結(jié)果集中，然后將與該節(jié)點存在連接的邊刪除掉，迭代至無向圖G中邊為空為止.

進(jìn)一步，實現(xiàn)了基于MapReduce的貪心策略GRDC.在Map處理階段，為了提高GRDC算法的運行效率，先將每行數(shù)據(jù)的對應(yīng)MACS集合進(jìn)行特征區(qū)域多樣性選擇，得到若干代表此行序列的MACS子集.因為當(dāng)特征維數(shù)比較高時，每行數(shù)據(jù)對應(yīng)求交得到的MACS數(shù)量相對也比較大，可以先在每行數(shù)據(jù)對應(yīng)的集合內(nèi)部進(jìn)行多樣性選擇，使得最后進(jìn)行多樣性選擇的MACS數(shù)量大大減少.然后，依次對原始數(shù)據(jù)集中每一行數(shù)據(jù)對應(yīng)的MACS集合都執(zhí)行此操作.在Reduce處理階段，先得到對應(yīng)Map任務(wù)處理后的結(jié)果，然后將匯總得到的MACS集合再進(jìn)行一次特征區(qū)域多樣性選擇，得到最終經(jīng)過多樣性選擇后的特征區(qū)域MACS集合，并將其寫入到HDFS中.

3.5 基于置換搜索的并行蟻群交互特征挖掘

本節(jié)將在分解后的特征區(qū)域內(nèi)，利用蟻群算法發(fā)現(xiàn)k階交互作用.

3.5.1 基本蟻群算法

本節(jié)引入面向k階交互的基本蟻群算法.其具體步驟如下：

首先，初始化antNum個螞蟻對象和迭代次數(shù)iterNum.然后，每個螞蟻隨機(jī)選擇一個初始特征，原始數(shù)據(jù)集中的每個特征對應(yīng)一個初始信息素值τ0.每一個螞蟻根據(jù)式(4)選擇下一個特征，直到選擇出的特征子集大小為k.比較每只螞蟻搜索得到的顯著特征子集，選擇出顯著性特征子集的并根據(jù)式(5)更新集合中每個特征的信息素濃度大小.

(4)

τk(t+1)=(1-ρ)τk(t)+Δτk(t),

(5)

3.5.2 置換策略

基本蟻群算法在每一次蟻群迭代過程中，得到的顯著交互特征子集信息不能夠保存下來.在隨后的迭代過程中，是根據(jù)每個特征位點的概率密度進(jìn)行重新選擇的.針對此問題，提出基于置換的蟻群算法，使得改進(jìn)后的蟻群優(yōu)化算法能夠避免顯著特征所帶來的邊際作用影響.基本思想為:保留上一步迭代過程中得到的顯著交互特征子集，然后每次對集合內(nèi)的一個特征進(jìn)行置換.在選擇下一個特征時，根據(jù)每個特征已置換次數(shù)選擇集合外的一個特征與集合內(nèi)的每個特征進(jìn)行置換，從而減少特征的重復(fù)選擇增加多樣性.如果新的特征子集較原來的顯著，就執(zhí)行此次置換，否則不執(zhí)行.這種方法，能夠很好地利用高階交互特征子集信息，有效防止低階特征所帶來的邊際作用影響，從而得到顯著的高階交互特征子集.本文中，該策略被記為RouteSearching，共包括3個階段：

1) 根據(jù)式(4)選擇大小為k的初始位點數(shù).在初始位點選擇過程中，根據(jù)計算每個位點對應(yīng)打分值Score，然后計算出對應(yīng)的顯著性p_value，最終選擇出大小為k的位點子集.

2) 對迭代次數(shù)以及是否執(zhí)行下次置換標(biāo)記進(jìn)行初始化，并通過置換搜索的思想，迭代地選出顯著位點子集，并將其全部加入到集合S′中.當(dāng)已選位點子集外沒有使該特征子集更顯著的位點時，迭代停止.

3) 對集合S′中的每個顯著位點子集進(jìn)行特征位點的后向檢測，目的是為了刪除掉冗余特征子集，最終得到顯著的不大于k階的交互特征子集集合S.

由于空間有限，此處略去RouteSearching的具體實現(xiàn)偽代碼.

3.5.3 基于置換的蟻群算法

在基本蟻群算法的迭代搜索中，加入置換策略RouteSearching后，得到改進(jìn)的蟻群算法，步驟如下：

步驟1. 計算數(shù)據(jù)集中每個位點的Score，選擇出顯著位點，將其移除. 這樣能夠減少強(qiáng)相關(guān)特征所帶來的邊際作用的影響.

步驟3. 通過置換策略RouteSearching選出顯著交互特征子集集合S′，并將其加入到結(jié)果集合S中.

步驟4. 對步驟3得到的集合S′中的每個顯著特征位點進(jìn)行信息素大小更新操作.然后進(jìn)入到下一個螞蟻迭代搜索中，直到所以螞蟻迭代搜索完成，算法終止.

3.5.4 并行蟻群算法

在本節(jié)中，基于MapReduce模型對基于置換的蟻群算法進(jìn)行了并行化，提出了相應(yīng)的算法PACOIFS.以下將給出詳細(xì)的描述.

蟻群優(yōu)化算法從本質(zhì)上來說就屬于一種并行算法，具有比較高的可并行性.在蟻群中的每只螞蟻搜索的過程其實就是得到可行解的過程，螞蟻之間搜索過程是相互獨立的，每只螞蟻搜索可行解的過程只與信息素濃度大小以及啟發(fā)式函數(shù)有關(guān)，并且當(dāng)整個種群迭代搜索完成之后，才進(jìn)行一次全局信息素濃度大小的更新操作.因此，可以把每個螞蟻放到不同的機(jī)器節(jié)點上并行地搜索可行解，這樣就使得集群并行化運行的優(yōu)勢得到了充分利用.蟻群算法比較常用的并行方式有2種：1)并行交互蟻群，即是將整個蟻群劃分為幾小部分，然后在這些部分蟻群之間并行；2)并行螞蟻，即是將整個蟻群中的每一個螞蟻進(jìn)行并行.通過以上介紹可以發(fā)現(xiàn)，并行螞蟻方式其實是并行交互蟻群方式的一種極端化形式，即是將整個蟻群劃分成螞蟻數(shù)目個小部分蟻群.由于并行交互蟻群方式每臺機(jī)器所對應(yīng)Mapper任務(wù)迭代時間會很長，在數(shù)據(jù)集規(guī)模比較大的情況下對機(jī)器性能要求較高.本文將采用并行螞蟻方式，因為這樣能夠?qū)⒂嬎闳蝿?wù)更好地分配到各個機(jī)器節(jié)點上，減輕每臺機(jī)器的計算壓力.偽碼如算法3所示.

算法3.PACOIFS().

Map Function(key,value=第k個螞蟻對應(yīng)的交互特征子集集合Sk).

輸入：SNP數(shù)據(jù)集、信息素以及蟻群參數(shù)等信息；

輸出：〈key,value=Sk〉.

① 讀取數(shù)據(jù)集D并劃分成m塊；

② for 每個塊bido

③ 讀取信息素及蟻群參數(shù)等信息進(jìn)行初始化;

④Sk←RouteSearching(bi,α,k,maxQ,SNPs);

⑤ 輸出〈key,value=Sk〉.

⑥ end for

Reduce Function(key,value=顯著交互特征子集集合S).

輸入:〈key,Sk〉；

輸出:〈key,S〉.

① for 在Map任務(wù)中的每個Skdo

② 添加SktoS;

③ end for

④ for 每個S中的SNP do

⑤ 根據(jù)式(5)更新每一個SNP位點的信息素;

⑥ end for

⑦ 將更新后的信息素向量信息寫入到HDFS之上;

⑧ 輸出〈key,S〉.

根據(jù)本文中所采用的蟻群并行方式可知，在MapReduce框架中，局部搜索的過程即對應(yīng)著Map接口，全局信息素濃度大小更新操作即對應(yīng)著Reduce接口.每個MapReduce Job即對應(yīng)一次蟻群進(jìn)行迭代搜索交互特征子集的過程，其中每個Mapper任務(wù)即對應(yīng)著一組螞蟻搜索交互特征子集的過程，有多少組螞蟻就會對應(yīng)多少個Mapper任務(wù)，每一個Reduce任務(wù)即是將Mapper任務(wù)結(jié)果匯總并寫入到HDFS中以及信息素更新操作的過程，整個蟻群的迭代即是MapReduce Job的迭代.由于在MapReduce框架下開啟Mapper或者Reduce會有相應(yīng)的自身開銷，本文將會使用少量的螞蟻數(shù)目，把整個算法的迭代搜索盡量放到每個Mapper任務(wù)中去，3.5.2節(jié)給出的置換策略恰恰符合這種需求，因為本文是基于每個螞蟻在迭代搜索過程中通過置換的思想來挖掘交互特征子集的.

4 實驗與分析

4.1 實驗環(huán)境

本實驗所用的軟硬件環(huán)境是集群由5臺服務(wù)器組成，每臺服務(wù)器的配置為Red hat 64位操作系統(tǒng)，16核CPU，主頻1.9 GHz，16 GB內(nèi)存，2 TB硬盤；Hadoop版本為2.6.0，Spark版本為1.6.0，Java版本為1.8.0，Scala版本為2.10.4.單機(jī)開發(fā)環(huán)境：操作系統(tǒng)為Windows7 32位旗艦版，主頻3.10 GHz，4 GB內(nèi)存，500 GB硬盤；開發(fā)工具IntelliJ IDEA Community Edition 15.0.2，Java版本為1.8.0，Scala版本為2.10.4.

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

1) 模擬數(shù)據(jù)集.采用由Marchini等人[19]提出的相加(模型1)、倍增(模型2)以及閾值(模型3)三種基因模型.通過設(shè)置各個模型中不同的參數(shù)值來產(chǎn)生相應(yīng)的實驗測試數(shù)據(jù).模擬數(shù)據(jù)集中共有4 000個樣本(其中包含2 000個病例和2 000個對照樣本)和2 000個SNP位點.由于數(shù)據(jù)的統(tǒng)計能力受MAF的影響較大，取其值為0.1，0.2，0.5；邊際作用λ取值較小，在模型1中設(shè)置大小為0.3，在模型2以及模型3中都設(shè)置為0.2；患病率(Prevalence)在每個模型中取值分別為0.1，0.4，0.7；對于連鎖不平衡r2大小，每個模型分別取值0.7和1.此外，還通過增加每個模型對應(yīng)數(shù)據(jù)集中的數(shù)據(jù)行數(shù)和SNP位點數(shù)來檢測算法在大數(shù)據(jù)集中的統(tǒng)計能力.

2) 真實數(shù)據(jù)集.使用了糖尿病患者數(shù)據(jù)作為真實數(shù)據(jù)集.該糖尿病患者數(shù)據(jù)集從美國凱斯西儲大學(xué)獲取.該數(shù)據(jù)集共包含有90個樣本個體，其中包括45個日本人和45個中國人.日本人中有38個患病(case)樣本，中國人中有11個患病(case)樣本，其他都為未患病(control)樣本，并且每個樣本都由22條染色體組成，每條染色體所對應(yīng)的SNP位點個數(shù)在3 000～20 000之間，詳細(xì)信息如表2所示：

Table 2 The Diabetes Dataset 表2 糖尿病數(shù)據(jù)集

Fig. 4 Power comparison on MAF change for different models圖4 不同模型的識別力Power受MAF變化比較圖

4.3 參數(shù)設(shè)置

本文算法中主要包含參數(shù)：統(tǒng)計顯著性閾值α、待選擇SNP子集大小k、局部搜索算法Route-Searching迭代次數(shù)上限值maxQ、蟻群種群大小antNum、種群迭代次數(shù)iterNum、SNP位點初始信息素大小τ0以及信息素?fù)]發(fā)系數(shù)ρ.

統(tǒng)計顯著性閾值α通過多重假設(shè)檢驗校正設(shè)置為α=0.01，每個SNP位點的初始信息素大小τ0=1，信息素?fù)]發(fā)系數(shù)ρ=0.05，待選擇SNP位點集合大小k∈[1,5].由于在置換搜索RouteSearching中循環(huán)條件由當(dāng)前迭代次數(shù)和是否執(zhí)行下次置換參數(shù)共同決定，為了防止每個螞蟻在局部置換搜索中迭代次數(shù)過多，并且又能夠利用已經(jīng)得到的顯著SNP位點集合信息，將最大迭代次數(shù)設(shè)置為maxQ=[0.01,0.2]，默認(rèn)值100.蟻群種群大小antNum∈[0.01,0.5]×N，種群迭代次數(shù)iterNum∈[0.05,0.7]×N，蟻群算法中2個重要的參數(shù)是種群大小antNum和種群迭代次數(shù)iterNum.

4.4 識別能力分析

本節(jié)將在模擬數(shù)據(jù)集上對提出的算法進(jìn)行識別能力分析.算法識別能力是指算法能夠正確檢測出致病位點的數(shù)據(jù)集個數(shù)在所有數(shù)據(jù)集中所占的比例大小，用式Power=NαN表示.其中，N表示每個模型中對應(yīng)的數(shù)據(jù)集的數(shù)目，Nα表示算法能夠正確識別出致病位點的數(shù)據(jù)集數(shù)目.根據(jù)4.3節(jié)中介紹的各基因模型參數(shù)設(shè)置，生成對應(yīng)的模擬數(shù)據(jù)集，每種模型對應(yīng)18種不同參數(shù)組合.實驗中將要對比的參數(shù)分別為次等位基因頻率MAF、患病率Prevalence以及連鎖不平衡r2.基于各模型對應(yīng)的數(shù)據(jù)集對文中算法PACOIFS進(jìn)行實驗測試，并將實驗結(jié)果和代表性算法SNPHarvester[13]、AntEpiSeeker[14]以及貝葉斯上位關(guān)聯(lián)映射(BEAM)[12]進(jìn)行對比分析.具體實驗結(jié)果分別見圖4～6所示.

圖4中，參數(shù)MAF的變化和模型的不同對各個算法的識別能力的影響還是較大的.其中，對于本文算法PACOIFS和AntEpiSeeker算法大致是隨著MAF的增大其Power隨之增大，表明當(dāng)增大致病SNP位點的次等位基因頻率，MAF會增大算法對致病SNP位點集合的檢測能力.這2個算法在3個模型中的識別能力大致相同，都可以達(dá)到不錯的效果.算法BEAM，SNPHarvester的波動在圖4(a)和圖4(c)下不太穩(wěn)定，在圖4(b)下是恰好呈相反的趨勢，即是隨著MAF的增大其Power反而減小，即是表明當(dāng)增大致病SNP位點的次等位基因頻率MAF會減小算法對致病SNP位點集合的檢測能力.

圖5中，患病率Prevalence的變化對算法Ant-EpiSeeker和BEAM的影響不太大，大致在圖5(b)和圖5(c)上呈輕微上升趨勢，而在圖5(a)中當(dāng)Prevalence=0.4時性能較差，總體上算法Ant-EpiSeeker要比算法BEAM的Power性能好一些.本文算法PACOIFS在各模型中則表現(xiàn)出和算法AntEpiSeeker相似的趨勢，其Power大致上比算法AntEpiSeeker好一些.算法SNPHarvester則表現(xiàn)出和其他算法不一樣的性能，其在圖5(a)和圖5(b)中當(dāng)Prevalence=0.4時性能較差，在圖5(c)中Power隨著Prevalence增大而呈上升趨勢.

Fig. 5 Power comparison on Prevalence change for different models圖5 不同模型的Power受患病率變化比較圖

圖6中，算法BEAM,SNPHarvester,AntEpi-Seeker以及本文算法PACOIFS在圖6(a)～(c)中的Power值都隨著連鎖不平衡r2的增大而提升.其中，算法AntEpiSeeker和本文算法PACOIFS表現(xiàn)較其他算法好，在圖6(b)中這2個算法都達(dá)到了最好值Power=1.在圖6(c)中的Power性能比圖6(a)和圖6(b)稍差，并且算法BEAM，SNPHarvester大致也呈這種趨勢.

Fig. 6 Power comparison on r2 change for different models圖6 不同模型的Power受r2變化比較圖

通過設(shè)置不同的參數(shù)，即蟻群種群大小和迭代次數(shù)，對本文算法PACOIFS進(jìn)行相應(yīng)的識別能力Power對比分析.首先，生成圖6(b)默認(rèn)參數(shù)下規(guī)模為2 000×2 000的模擬數(shù)據(jù)集，然后設(shè)置種群大小antNum=100，iterNum分別為200，250，300，350，400，迭代次數(shù)iterNum=200，antNum分別設(shè)置為100，150，200，250，300，將實驗結(jié)果和算法AntEpiSeeker在相同實驗參數(shù)下進(jìn)行對比分析.結(jié)果如圖7～8所示：

Fig. 7 Power comparison on iterNum change圖7 Power受iterNum變化比較圖

Fig. 8 Power comparison on antNum change圖8 Power受antNum變化比較圖

4.5 可擴(kuò)展性分析

本節(jié)將利用模擬數(shù)據(jù)集對算法PACOIFS進(jìn)行可擴(kuò)展性方面的分析.首先，觀察分析提出算法隨數(shù)據(jù)行數(shù)及SNP位點數(shù)的增加，其Power性能的變化，并將實驗結(jié)果與SNPHarvester[13]，AntEpiSeeker[14]，BEAM[12]進(jìn)行對比.基于基因模型2的實驗結(jié)果分別如圖9和圖10所示.

Fig. 9 Power comparison on datalines change圖9 Power受數(shù)據(jù)行數(shù)變化比較圖

Fig. 10 Power comparison on SNP-locus change圖10 Power受SNP位點數(shù)變化比較圖

圖9中，算法BEAM,SNPHarvester,AntEpi-Seeker及算法PACOIFS在模型2中的Power性能隨數(shù)據(jù)行數(shù)的增加而增大.其中，算法AntEpi-Seeker和本文算法PACOIFS在SNP位點數(shù)為2 000、數(shù)據(jù)行數(shù)為3 000和4 000時均達(dá)到了最好性能.因為隨著樣本數(shù)的增多其統(tǒng)計能力也會隨著增大，并且隨著數(shù)據(jù)行數(shù)的增加本文算法要比算法AntEpiSeeker的總體Power性能要好.圖10中，算法BEAM,SNPHarvester,AntEpiSeeker以及算法PACOIFS在模型2中的Power性能隨著SNP位點數(shù)增加而大致減小，其中算法AntEpiSeeker和算法PACOIFS在數(shù)據(jù)行數(shù)固定為2 000對應(yīng)SNP位點數(shù)為100的時候均達(dá)到了最好性能.當(dāng)SNP位點數(shù)逐漸增加的時候本文算法PACOIFS也是呈最好的Power性能，因為本文在算法框架加入了特征區(qū)域劃分階段，使得SNP位點的增加對其算法性能影響不大.這也是大多數(shù)算法在高維數(shù)據(jù)中性能表現(xiàn)不好的原因.

算法隨著SNP位點數(shù)增加的響應(yīng)時間變化情況如圖11所示.由于真實情況下，數(shù)據(jù)集中的樣本量一般很少而SNP位點數(shù)一般會很大，所以這里只分析隨著SNP位點數(shù)的增加，各算法在時間效率上的變化情況.同樣是在基因模型2上采用不同位點數(shù)進(jìn)行實際測試，其中樣本數(shù)為2 000.由圖11可見，算法BEAM和SNPHarvester相對比較快，而算法AntEpiSeeker和本文算法較慢.這是因為算法AntEpiSeeker加入了后續(xù)最小化假陽性處理階段，并且其使用了2階段蟻群迭代，所以時間相對較長.但是，本文算法ACOIFS由于加入了無關(guān)特征削減、序列之間求MACS等前期處理過程，所以運行時間也比較長，這也是其Power性能比較好的原因.隨著SNP位點數(shù)增加，算法BEAM和SNPHarvester的響應(yīng)時間增長幅度較明顯，而算法AntEpiSeeker和本文算法的響應(yīng)時間雖然也在增加，但并不明顯.并行算法PACOIFS表現(xiàn)最好，在SNP位點數(shù)為1 000的時候其運行時間已經(jīng)達(dá)到最少，并且隨著SNP位點數(shù)增加，這種趨勢更加明顯.

Fig. 11 Runtime comparison on SNP-locus change圖11 運行時間受SNP位點數(shù)變化比較圖

4.6 加速比分析

加速比是同一個任務(wù)在不同處理器或并行處理器中運行時間的比值，用來衡量程序代碼并行化或并行系統(tǒng)的真實有效性.在獨占CPU資源的前提下，假設(shè)某個串行程序運行時間為Ts，而將該程序并行化后，k個進(jìn)程在k個CPU并行運行的時間為Tk，那么程序并行化之后的加速比R則為R=TsTk.由于在本文算法框架中包括特征區(qū)域劃分以及特征區(qū)域選擇階段，為的是將原始長的序列切分成若干短的序列并對其進(jìn)行多樣性選擇，這樣同時也使得數(shù)據(jù)易于并行化處理.將在真實數(shù)據(jù)集上采用不同SNP位點數(shù)來測試算法PACOIFS基于MapReduce框架的加速比，結(jié)果見圖12所示:

Fig. 12 Speed up on different datasets and nodes圖12 不同數(shù)據(jù)量不同節(jié)點個數(shù)的加速比

由圖12可見，本文算法的加速比隨著節(jié)點個數(shù)的增加逐漸增加，同時增加的趨勢在不斷減小.也就是說，如果不斷地增加節(jié)點，則處理時間不會呈線性的減少.這主要是因為MapReduce框架下，數(shù)據(jù)傳輸會消耗時間.由MapReduce模型本身的特點可知，Map任務(wù)完成后，需要將處理結(jié)果進(jìn)行排序，并根據(jù)partition函數(shù)將數(shù)據(jù)傳輸?shù)綄?yīng)的Reduce節(jié)點上，即shuffle處理階段，這段時間需要進(jìn)行數(shù)據(jù)拷貝.當(dāng)數(shù)據(jù)集較大時，增加節(jié)點雖然能使數(shù)據(jù)處理的時間減少，但并沒使得數(shù)據(jù)傳輸?shù)臅r間減少.因此，加速比曲線沒有呈斜率為1的趨勢增加.

4.7 有效性分析

根據(jù)以上實驗結(jié)果，算法AntEpiSeeker和本文算法PACOIFS性能表現(xiàn)最接近，因為AntEpiSeeker使用了2階段蟻群迭代搜索并加入了后續(xù)最小化假陽性處理階段，而本文算法由于加入了無關(guān)特征削減以及特征區(qū)域劃分求MACS等前期處理過程，所以實驗結(jié)果相對較好.在算法有效性分析中，將基于糖尿病患者數(shù)據(jù)集對PACOIFS和AntEpiSeeker進(jìn)行比較測試和對比分析.首先，根據(jù)算法AntEpi-Seeker中的建議將實驗參數(shù)設(shè)置為iItCountLarge=10 000，iItCountSmall=20 000，其他參數(shù)取默認(rèn)值；本文算法參數(shù)設(shè)置為antNum=2 000，iterNum=20 000.

表3中列舉出了2種算法識別出的部分顯著位點集合.可以看出，本文算法PACOIFS所識別出的顯著位點集合比算法AntEpiSeeker數(shù)目多，而且結(jié)果更顯著.由于本文算法在開始時將比較顯著的單SNP位點提前移除掉，因此得到的結(jié)果更好地反映出了在全基因組范圍內(nèi)SNP位點交互作用分布的情況.最后，將選擇出的結(jié)果和糖尿病的參照致病位點集合進(jìn)行對比，進(jìn)一步驗證了結(jié)果的可靠性.

Table 3 Some Significant SNP Groups Identified by AntEpiSeeker and PACOIFS表3 AntEpiSeeker和PACOIFS識別出的部分顯著SNP位點集合

5 結(jié) 論

針對大規(guī)模生物序列中高階交互特征挖掘面臨的密集計算問題，本文提出了一種新的基于并行處理和演化計算的解決框架.與現(xiàn)有框架不同，提出的框架從減少分組數(shù)據(jù)位點數(shù)入手，對數(shù)據(jù)進(jìn)行垂直分解.由于位點數(shù)是制約互作挖掘效率的根本因素，提出的框架具有天然的效率優(yōu)勢.同時，提出了極大等位公共子序列MACS的概念，能保證劃分后交互作用的“高內(nèi)聚低耦合”性，避免了分解后計算任務(wù)的耦合，減少了結(jié)果的丟失.在模擬數(shù)據(jù)集和真實的糖尿病患者數(shù)據(jù)集上，從識別能力、可擴(kuò)展性、加速比、有效性等角度，將PACOIFS與相關(guān)經(jīng)典算法進(jìn)行了比較分析，結(jié)果表明提出的算法具有明顯的效率和有效性優(yōu)勢.由于本文提出的算法框架采用求MACS的方式對特征區(qū)域進(jìn)行劃分操作，當(dāng)樣本較多時比較耗時，今后可以考慮利用已得到的結(jié)果減少求MACS的操作.另一方面，本文是在得到所有求MACS后進(jìn)行多樣性特征區(qū)域的選擇，未來可以考慮選擇一種在線處理方式進(jìn)行多樣性特征區(qū)域選擇，以提高效率.