王 超 英

(東莞職業(yè)技術(shù)學(xué)院 廣東 東莞 523808)

0 引 言

隨著醫(yī)院的信息化建設(shè)和普及，全國(guó)的醫(yī)療系統(tǒng)每天產(chǎn)生海量的數(shù)據(jù)，如何高效地管理醫(yī)療大數(shù)據(jù)成為了一個(gè)難題[1]。在醫(yī)學(xué)領(lǐng)域存在大量的高維小樣本數(shù)據(jù)，例如醫(yī)囑類(lèi)文檔數(shù)據(jù)、病例類(lèi)文檔數(shù)據(jù)、DNA微陣列數(shù)據(jù)[2-3]。海量高維小樣本數(shù)據(jù)集的聚類(lèi)和分析需要大量的處理時(shí)間。以DNA微陣列數(shù)據(jù)為例，微陣列的維度高，且含有大量的冗余樣本和不相關(guān)樣本，嚴(yán)重影響聚類(lèi)器的準(zhǔn)確率[4]。

通過(guò)特征選擇處理能夠降低特征集的維度，刪除冗余特征和不相關(guān)特征，有助于提高數(shù)據(jù)聚類(lèi)的準(zhǔn)確率和速度，是當(dāng)前高維小樣本數(shù)據(jù)挖掘領(lǐng)域的研究熱點(diǎn)[5]。為了解決上述高維小樣本數(shù)據(jù)集的瓶頸問(wèn)題，許多研究人員進(jìn)行了深入的分析[6]。文獻(xiàn)[7]基于隨機(jī)森林變量重要性分?jǐn)?shù)提出了一種加權(quán)K-均值的基因微陣列數(shù)據(jù)聚類(lèi)算法，以基因?yàn)閷?duì)象、樣本為特征、基因的重要性分?jǐn)?shù)為權(quán)重進(jìn)行K-均值聚類(lèi)，該算法提高了聚類(lèi)結(jié)果的同質(zhì)性和差異性。文獻(xiàn)[8]提出的高維微陣列數(shù)據(jù)混合特征選擇算法分為兩層：第一層使用信噪比方法計(jì)算全部基因的信噪比值，根據(jù)信噪比值選擇指定數(shù)目的信息基因，過(guò)濾無(wú)關(guān)基因；第二層使用改進(jìn)的Lasso方法對(duì)第一層得到的信息基因候選子集進(jìn)行特征選擇，剔除冗余基因。文獻(xiàn)[7-8]均采用了過(guò)濾式特征選擇算法，分別將隨機(jī)森林變量重要性分?jǐn)?shù)和基因信噪比作為特征重要性的評(píng)價(jià)指標(biāo)，此類(lèi)傳統(tǒng)方法難以獲得最少的特征數(shù)量[9]。

DNA微陣列數(shù)據(jù)存在以下特點(diǎn)[10]：(1) 生物標(biāo)志物的數(shù)量極低(低于1%);(2) 生物標(biāo)志物能夠提高特征子集的聚類(lèi)準(zhǔn)確率;(3) 過(guò)濾式特征選擇可縮小搜索空間，但無(wú)法保證優(yōu)質(zhì)的特征選擇結(jié)果;(4) 樣本量少，聚類(lèi)器訓(xùn)練難度大。提取最少量、最高準(zhǔn)確率的特征子集是特征選擇的兩個(gè)目標(biāo)，采用元啟發(fā)式算法對(duì)兩個(gè)目標(biāo)進(jìn)行尋優(yōu)是一個(gè)重要的研究方向，例如：量子進(jìn)化算法[11]、飛蛾撲火優(yōu)化算法[12]和人工蜂群優(yōu)化算法[13]，但此類(lèi)算法在全局搜索和局部開(kāi)發(fā)之間未能實(shí)現(xiàn)較好的權(quán)衡，同時(shí)基于迭代的算法對(duì)于海量數(shù)據(jù)的處理速度較慢，效率較低。文化基因算法[14](Memetic Algorithm,MA)是一種基于種群全局搜索和基于個(gè)體局部啟發(fā)式開(kāi)發(fā)的結(jié)合體，其搜索效率最高比遺傳算法快數(shù)個(gè)數(shù)量級(jí)。鑒于遞歸模型[15]在全局搜索和局部開(kāi)發(fā)之間能夠?qū)崿F(xiàn)較好的平衡，本文將文化基因算法和遞歸模型結(jié)合，設(shè)計(jì)了遞歸文化基因算法(Recursive Memetic Algorithm,RMA)。采用遞歸文化基因算法求解高維小樣本數(shù)據(jù)的特征選擇問(wèn)題，但是基于迭代的尋優(yōu)過(guò)程依然較為耗時(shí)。

分布式并行計(jì)算模型是加快大數(shù)據(jù)處理的有效方法，目前的大數(shù)據(jù)技術(shù)大多基于MapReduce框架，但該框架并不適合基于迭代的算法，原因主要有[16]：(1) 每次迭代作為獨(dú)立任務(wù)重新處理，需要重新讀寫(xiě)數(shù)據(jù)和初始化數(shù)據(jù)；(2) 每次迭代需要重新從分布式文件系統(tǒng)加載程序，消耗I/O資源和CPU資源；(3) 前一次迭代必須全部結(jié)束并且數(shù)據(jù)全部寫(xiě)入分布式文件系統(tǒng)，才能開(kāi)始下一次迭代。Spark框架[17]則是基于內(nèi)存的分布式并行計(jì)算模型，將數(shù)據(jù)緩存于各節(jié)點(diǎn)的內(nèi)存中，從而縮短數(shù)據(jù)訪(fǎng)問(wèn)的延遲。本文基于Spark架構(gòu)設(shè)計(jì)和實(shí)現(xiàn)了遞歸文化基因算法，可以高效、準(zhǔn)確地進(jìn)行高維小樣本數(shù)據(jù)集的特征選擇處理。并且本文研究并提出了一種基于猶豫模糊集理論(Hesitant Fuzzy Set,HFS)的加權(quán)相關(guān)性指數(shù)度量基因的相似性，設(shè)計(jì)了基于HFS的密度聚類(lèi)方法，HFS對(duì)噪聲具有魯棒性，設(shè)計(jì)了迭代Spark框架的彈性分布式數(shù)據(jù)集(Resilient Distributed Datasets,RDD)模塊將基因集聚類(lèi)處理。

1 猶豫模糊集理論

HFS是一種擴(kuò)展模糊集，其隸屬度不是確定值或服從確定分布，而是多個(gè)可能值，該理論對(duì)于群決策信息的處理具有優(yōu)勢(shì)[18]。

1.1 HFS定義

設(shè)X={x1,x2,…,xn}為一個(gè)固定集，如果A滿(mǎn)足以下條件，則認(rèn)為A是X的一個(gè)HFS。

A={〈x,hA(x)〉|x∈X}

(1)

設(shè)h、h1、h2為3個(gè)HFE，HFE支持以下三個(gè)運(yùn)算：

(2)

(3)

(4)

給定h∈HFE(x)，x的下限和上限分別定義為：

LB:h-(x)=min(h(x))UB:h+(x)=max(h(x))

(5)

設(shè)Aenv(h)表示h的包絡(luò)線(xiàn)，表示為(h-,1-h+)。

1.2 HFS的相關(guān)性

設(shè)X={x1,x2,…,xn}為一個(gè)離散論域，A和B是X上的兩個(gè)HFS，表示為A={〈xi,hA(xi)〉|xi∈X,i=1,2,…,n}，B={〈xi,hB(xi)〉|xi∈X,i=1,2,…,n}。如果len(hA(xi))

(6)

兩個(gè)HFS的相關(guān)性定義為：

(7)

如果A,B∈HFS，那么CHFS(A,A)=EHFS(A),CHFS(A,B)=CHFS(B,A)。最終可計(jì)算出兩個(gè)HFS的相關(guān)系數(shù)：

(8)

HFS的相關(guān)系數(shù)具有以下屬性[19]：ρHFS(A,B)=ρHFS(B,A)，0≤ρHFS(A,B)≤1，如果A=B，那么ρHFS(A,B)=1。

2 聚類(lèi)系統(tǒng)總體結(jié)構(gòu)

本文的高維大數(shù)據(jù)聚類(lèi)系統(tǒng)(見(jiàn)圖1)主要分為兩個(gè)部分：基于密度的聚類(lèi)處理和遞歸文化基因的特征歸簡(jiǎn)處理。兩個(gè)部分是遞歸迭代的關(guān)系。在每次迭代中，遞歸文化基因算法遞歸地化簡(jiǎn)特征子集，將特征子集輸入聚類(lèi)處理部分，根據(jù)聚類(lèi)結(jié)果判斷特征子集的質(zhì)量，實(shí)現(xiàn)了一種封裝式特征選擇方式。密度聚類(lèi)處理應(yīng)用特征子集做聚類(lèi)處理，基于Spark框架提高處理的速度和效率。

圖1 高維大數(shù)據(jù)聚類(lèi)系統(tǒng)的主要模塊圖

3 基于遞歸文化基因的特征歸簡(jiǎn)

遞歸文化基因算法(RMA)進(jìn)行特征選擇總體上重復(fù)執(zhí)行文化基因算法(MA)，逐步縮小特征空間。MA采用精英機(jī)制，隨機(jī)建立種群，按聚類(lèi)的適應(yīng)度值(準(zhǔn)確率)將染色體排列。首先，每個(gè)染色體做局部搜索處理，將其替換為更好的染色體；然后，采用輪盤(pán)賭選擇交叉算子的染色體，應(yīng)用MA的交叉算子和變異算子。MA算法對(duì)于一般數(shù)據(jù)的特征選擇效果較好，但對(duì)高維微陣列數(shù)據(jù)的降維能力較弱，因此設(shè)計(jì)了RMA模型。

應(yīng)用MA優(yōu)化種群，選出適應(yīng)度高于動(dòng)態(tài)閾值θ的i個(gè)最優(yōu)染色體，i個(gè)染色體和β%的最優(yōu)特征合并組成新的特征空間，新特征子集的聚類(lèi)準(zhǔn)確率應(yīng)當(dāng)高于θ;再次應(yīng)用MA處理上述特征子集，逐漸縮小特征空間，并且提高后續(xù)基于密度聚類(lèi)的準(zhǔn)確率。圖2為RMA的算法流程。

圖2 RMA的算法流程

(1) 局部搜索。初始化階段將排序的特征集輸入RMA做局部搜索處理，局部搜索包括ADD和DEL兩個(gè)算子。設(shè)cr為選擇的一個(gè)基因子集的染色體，P和Q分別為染色體cr選擇和排除的基因子集。ADD算子將Q的基因加入P中，DEL算子將基因從P放入Q中。局部搜索的主要問(wèn)題是決定ADD、DEL算子的輸入基因集，首先采用ReliefF方法[20]將所有的特征排序，然后ADD算子選擇Q中最優(yōu)的特征，加入P中，DEL算子選擇P中最差的特征，移至Q中。

(2) 交叉算子。交叉算子是一種遺傳算子，交叉算子的方式主要有單點(diǎn)交叉、兩點(diǎn)交叉、多點(diǎn)交叉、融合交叉和均勻交叉等，其中兩點(diǎn)交叉實(shí)現(xiàn)較為簡(jiǎn)單，且無(wú)須設(shè)置交叉概率，因此采用兩點(diǎn)交叉算子。

(3) 變異算子。變異算子也是一種遺傳算子，有助于增強(qiáng)算法的全局搜索能力。變異率設(shè)為隨機(jī)變量q，每次迭代基于概率q應(yīng)用變異算子。

(4) 綜合多目標(biāo)的適應(yīng)度函數(shù)?；蛭㈥嚵械木垲?lèi)準(zhǔn)確率結(jié)果作為主目標(biāo)，特征數(shù)量作為次目標(biāo)。具體方法為：設(shè)置一個(gè)誤差域值，在最大化次目標(biāo)的情況下，保留其中準(zhǔn)確率誤差小于閾值的染色體，刪除高于閾值的染色體。算法1所示為適應(yīng)度評(píng)價(jià)算法，采用加權(quán)平均準(zhǔn)確率和特征數(shù)量決定優(yōu)質(zhì)的染色體，其中γ=1-(選擇的特征數(shù)量/特征總數(shù)量)。

算法1適應(yīng)度評(píng)價(jià)算法

輸入：染色體ci,cj。

輸出：優(yōu)質(zhì)染色體c。

1. if ((mod(acc(ci)-acc(cj))>α))

//mod表示模運(yùn)算，acc()為準(zhǔn)確率

2. if (acc(ci) >acc(cj))

//比較準(zhǔn)確率

3. returnci;

4. else

5. returncj;

6. else

7.v=((w1×acc(ci))+(w2×γ(ci)))-((w1×acc(cj))+(w2×γ(cj)))

8. if (v>0)

9. returnci;

10. else

11. returncj

(5) 特征空間化簡(jiǎn)。設(shè)一個(gè)動(dòng)態(tài)閾值φ，種群的top-i染色體與φ作比較，如果染色體的適應(yīng)度高于φ，則重建種群。通過(guò)調(diào)節(jié)φ將種群規(guī)模從100%調(diào)節(jié)至ρ%，重建種群時(shí)，φ設(shè)為top-i種群的平均值。算法2是φ的增加算法偽代碼，其中：r是種群重建的次數(shù)，φ隨著r遞增，參數(shù)l1,l2,l3,step等參數(shù)根據(jù)實(shí)驗(yàn)確定。算法3是φ的衰減算法偽代碼，φ隨著recn增加而降低。如果處理測(cè)試集的特征集較小且準(zhǔn)確率高，或者達(dá)到了預(yù)定義的最大遞歸次數(shù)，則停止遞歸模型。根據(jù)實(shí)驗(yàn)結(jié)果將li值設(shè)為：l1=8,l2=6,l3=3,l4=10,l5=13。

算法2φ的遞增算法

輸入：top-i染色體的平均準(zhǔn)確率：accavg。

輸出：動(dòng)態(tài)閾值φ。

1.stepl=step×2；

2. if (r>l1)

3.φ=max(accavg,100-stepl);

4. else if (r>l2)

5.φ=max(accavg,100-2×stepl);

6. else if (r>l3)

7.φ=max(accavg,100-3×stepl);

8. end if

算法3φ的衰減算法

輸入：種群重建次數(shù)：recn。

輸出：動(dòng)態(tài)閾值φ。

試驗(yàn)過(guò)程中，每周從每個(gè)重復(fù)組中隨機(jī)選出1只小鼠，宰殺解剖后分別切取十二指腸段、空腸段、回腸段各3 cm左右。按照ELISA試劑盒使用說(shuō)明書(shū)通過(guò)酶標(biāo)儀進(jìn)行測(cè)定腸道黏膜SIgA的表達(dá)水平。試劑盒購(gòu)置于德國(guó)IBI分裝公司。

1. if (recn>l4&&φ≥ρ+step)

2.φ=φ-step;

3. if (recn>l5&&φ≥ρ-step)

4.φ=φ-step;

5. end if

4 基于密度的聚類(lèi)系統(tǒng)

4.1 聚類(lèi)算法

設(shè)ξ為類(lèi)的最少成員數(shù)量，τ為類(lèi)內(nèi)成員之間的最大距離。本文的研究對(duì)象為DNA微陣列數(shù)據(jù)集，所以假設(shè)數(shù)據(jù)元素為基因數(shù)據(jù)。設(shè)Nτ(gi)表示與基因(gi)距離小于τ的基因集，如果(|Nτ(gi)|+1)≥ξ，則稱(chēng)gi為一個(gè)核心基因(Core Gene,CG)，|Nτ(gi)|表示Nτ(gi)的大小。類(lèi)中的非CG基因稱(chēng)為邊緣基因(Border Gene,BG)。設(shè)G={g1,g2,…,gn}為一個(gè)基因集，如果gb∈CG且ga∈Nτ(gb)，那么ga由gb密度直達(dá)，將gb由ga密度直達(dá)表示為ga→gb，核心基因具有自我密度直達(dá)的性質(zhì)。設(shè)一個(gè)基因序列為{g1,g2,…,gk}，如果存在一個(gè)密度直達(dá)序列滿(mǎn)足條件gi→gi+1，1≤i≤k-1，那么g1向gk密度可達(dá)，表示為g1?gk。密度可達(dá)具有對(duì)稱(chēng)性，即g1?gk≡gk?g1。如果一個(gè)非CG與任意的CG均不是密度可達(dá)，那么該基因?yàn)楣铝⒒?Outlier Gene,OG)。CG的鄰居基因之間密度相連，并且所有可達(dá)基因均與CG密度相連。如果存在gk?gi，gk?gj，gk∈CG，那么gi和gj密度相連，本文將密度相連表示為gi?gj。對(duì)于聚類(lèi)問(wèn)題，同一個(gè)類(lèi)的基因之間密度相連。

4.2 基于猶豫模糊理論的基因相似性度量

相似性度量是決定聚類(lèi)準(zhǔn)確性的關(guān)鍵因素，采用HFS相關(guān)系數(shù)度量基因的相似性。設(shè)Gj為m個(gè)HFS的元素，cor=(ρij)m×m為一個(gè)相關(guān)矩陣。根據(jù)式(8)，ρij表示Gi和Gj間的相關(guān)系數(shù)，ρij具有以下屬性[19]：(1) 0≤ρij≤1,i,j=1,2,…,m；(2)ρii=1,i=1,2,…,m；(3)ρij=ρji,i,j=1,2,…,m。

Gk={〈Col,hGk(Col)〉|Col∈Co,Gk∈G}

(9)

將表達(dá)譜值Exp做歸一化處理，即{0≤Exp≤1|Exp?hGk(Col)}，hGk(Col)是基因Gk在條件hGk(Col)下的猶豫模糊元素，表示為HFE(k,l)。

(10)

WCHFS(Ga,Gb)=

(11)

(12)

5 分布式聚類(lèi)系統(tǒng)設(shè)計(jì)

采用基于密度的聚類(lèi)算法對(duì)基因微陣列做聚類(lèi)處理，基于Spark框架實(shí)現(xiàn)對(duì)聚類(lèi)算法的分布式處理?；赟park框架的分布式聚類(lèi)算法如圖1所示，主要分為5步：① 數(shù)據(jù)預(yù)處理；② 特征歸簡(jiǎn)；③ 計(jì)算距離；④ 建立子類(lèi)；⑤ 最終融合。

Spark框架由1個(gè)master模塊和若干的worker模塊組成。Master為worker分配job，在driver的協(xié)助下控制計(jì)算程序。worker的job從RDD、分布式文件系統(tǒng)(Hadoop Distributed File System,HDFS)、數(shù)據(jù)幀(DataFrame)等存儲(chǔ)系統(tǒng)讀取數(shù)據(jù)，然后處理剩余數(shù)據(jù)并輸出結(jié)果。worker節(jié)點(diǎn)的每個(gè)job由多個(gè)階段組成，worker串行地完成各個(gè)階段，每個(gè)階段可能獨(dú)立或者依賴(lài)之前階段的輸出。Spark內(nèi)部將數(shù)據(jù)分割為若干的分區(qū)，每個(gè)worker節(jié)點(diǎn)處理有權(quán)限的部分分區(qū)，輸出相應(yīng)的結(jié)果。worker節(jié)點(diǎn)度量數(shù)據(jù)的相似性，Spark的平臺(tái)負(fù)責(zé)實(shí)現(xiàn)分布式處理程序的細(xì)節(jié)，例如：負(fù)載均衡、容錯(cuò)機(jī)制、數(shù)據(jù)分布、尋址處理等。

Spark可獨(dú)立運(yùn)行或者建立于Hadoop YARN、Apache Mesos和Amazon EC2等云計(jì)算平臺(tái)。Spark支持分布式數(shù)據(jù)存儲(chǔ)系統(tǒng)，包括HDFS、HBase等。本文在YARN和EC2上建立Spark，采用HDFS作為分布式數(shù)據(jù)存儲(chǔ)系統(tǒng)。

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

基因表達(dá)數(shù)據(jù)集的預(yù)處理內(nèi)容主要為：采用表達(dá)譜對(duì)數(shù)變換處理[21]將輸入數(shù)據(jù)集歸一化，將一些數(shù)據(jù)做單位轉(zhuǎn)換和偏差轉(zhuǎn)換處理。如果實(shí)驗(yàn)矩陣40%以上的行是空值，那么過(guò)濾這種表達(dá)譜矩陣。例如：如果{(Ga,Gb)∈HFS,leni(hGa(xi))

每個(gè)微陣列實(shí)驗(yàn)是一個(gè)m×n的矩陣M，m是基因的數(shù)量，G=(G1,G2,…,Gm)，n是實(shí)驗(yàn)條件的數(shù)量，Co=(Co1,Co2,…,Con)。假設(shè)共存在q個(gè)醫(yī)學(xué)實(shí)驗(yàn)，最終共有n個(gè)條件下對(duì)m個(gè)基因的q次實(shí)驗(yàn)結(jié)果。設(shè)置一個(gè)Spark driver程序歸一化矩陣每行的數(shù)據(jù)，計(jì)算每行數(shù)據(jù)的均值μ和方差σ。將舊的表達(dá)譜替換為新值，計(jì)算式表示為：

(13)

5.2 特征歸簡(jiǎn)

采用第3節(jié)中基于遞歸文化基因的特征歸簡(jiǎn)算法縮小特征空間。

5.3 計(jì)算距離

加權(quán)的HFS相關(guān)系數(shù)WρHFS(Ga,Gb)表示基因Ga和Gb之間行為的相似性，相關(guān)系數(shù)越高說(shuō)明基因行為的相似度越高。從HDFS讀取基因表達(dá)譜矩陣，然后轉(zhuǎn)化為RDD格式。每個(gè)worker節(jié)點(diǎn)運(yùn)行一個(gè)迭代程序，并行地計(jì)算全部基因的加權(quán)信息能量WEHFS(Gk)。最終，通過(guò)key函數(shù)的reduce累加全部的能量值，將基因的ID作為key，信息量是key操作的目標(biāo)。計(jì)算加權(quán)相關(guān)系數(shù)之前，driver生成一個(gè)三角矩陣，稱(chēng)為watch列表。master節(jié)點(diǎn)檢查watch列表，驗(yàn)證兩個(gè)條件：(1) 計(jì)算不同基因集{Ga,Gb}之間的加權(quán)相關(guān)性；(2) 每一對(duì)基因的加權(quán)相關(guān)性應(yīng)當(dāng)僅計(jì)算一次，即WCHFS(Ga,Gb)=WCHFS(Gb,Ga)。

算法4所示是計(jì)算距離的偽代碼。在初始化階段為watch列表的每個(gè)元素設(shè)置一個(gè)預(yù)定值，作為是否已經(jīng)計(jì)算的標(biāo)記。逐漸將watch列表的元素替換為計(jì)算的加權(quán)相關(guān)系數(shù)，每個(gè)worker節(jié)點(diǎn)并行地計(jì)算基因的加權(quán)相關(guān)系數(shù)。最終master節(jié)點(diǎn)在每次迭代結(jié)束更新watch列表，master節(jié)點(diǎn)在watch列表的約束下修改RDD的內(nèi)容。在watch列表中沒(méi)有相關(guān)系數(shù)記錄的基因?qū)Ρ４嬗谙嗤腞DD中，然后計(jì)算它們的加權(quán)相關(guān)性并保存于watch列表中。重復(fù)該迭代程序，最終使所有基因?qū)哂邢嚓P(guān)性記錄。該程序具有可擴(kuò)展性，能夠處理任意規(guī)模的輸入數(shù)據(jù)矩陣。

算法4計(jì)算距離的偽代碼

輸入：watch列表WLE。

1. for each (WLEi?WLE)

2. WLEi←0;

//初始化為缺省值

3. end for

4. while ((WLEi?WLE)==0)

5. shuffle_rdd()；

6. 計(jì)算WCHFS；

7. 更新WLEi；

8. 計(jì)算WρHFS；

9. 生成GPM；

10. 建立DRM

算法5是建立密度可達(dá)矩陣的偽代碼。由近似矩陣生成密度可達(dá)矩陣，因?yàn)槊芏瓤蛇_(dá)具有對(duì)稱(chēng)性，所以無(wú)須保存對(duì)稱(chēng)的可達(dá)性，密度可達(dá)也是三角矩陣形式?；蚍植加诟鱾€(gè)worker節(jié)點(diǎn)，近似矩陣中基因Gi的第i列和第m+1-i行元素加入一個(gè)向量中，稱(chēng)為Vi，每個(gè)向量表示一個(gè)基因和其他所有基因的相似性，該向量被分配至各個(gè)worker節(jié)點(diǎn)。設(shè)計(jì)一個(gè)函數(shù)記錄Vi中加權(quán)相關(guān)系數(shù)高于閾值((WρHFS)i,j>τ)的基因數(shù)量Gj。每個(gè)節(jié)點(diǎn)并行地調(diào)用計(jì)數(shù)函數(shù)，將計(jì)數(shù)值高于ξ(|Nτ(Gi)|>ξ)的基因標(biāo)記為核心基因。

算法5建立密度可達(dá)矩陣

1. for each (Gi?RDD)

2.CntGi=0；

//初始化計(jì)數(shù)器

3. for each (Gj∈vector(Gi),i≠j)

4. if [(WρHFS)i,j>τ]

5.CntGi=CntGi+1；

6. end for

7. end for

8. for each (Gi?RDD)

9. if (CntGi>ξ)

10.Gi.coregene=TRUE；

//Gi標(biāo)記為核心基因

11. end for

12. fore ach (Gi?RDD &&Gi.coregene == TRUE)

13. for each(Gj∈vector(Gi),i≠j)

14. if ((WρHFS)i,j>τ)

15.Gi?Gj；

//更新RGL(Gi),RGL(Gj)

16. end for

17. end for

18. for each (Gi?RDD)

19. for each (Gj∈RGL(Gi),i≠j)

20. list(Gi)?list(Gj)；

//更新RGL()

21. end for

22. end for

標(biāo)記完所有的核心基因之后，將所有加權(quán)相關(guān)系高于閾值τ的基因標(biāo)記為對(duì)Gi直接密度可達(dá)，即{Gi→Gj|(Gj)∈Nτ(Gi),Gi∈CG}?？赡苣硞€(gè)基因?qū)τ诙鄠€(gè)基因密度可達(dá)，所以采用一個(gè)可達(dá)基因列表(Reachable Genes List,RGL)記錄對(duì)某個(gè)基因的密度可達(dá)基因。Gi被保存于Gj的可達(dá)基因列表中，Gj被保存于Gi的可達(dá)基因列表中。計(jì)算所有核心基因的直接密度可達(dá)之后，將直接密度可達(dá)泛化為密度可達(dá)。如果兩個(gè)基因之間密度直達(dá)，那么列表內(nèi)所有基因均為間接密度可達(dá)。

5.4 子聚類(lèi)階段

為每個(gè)基因Gi創(chuàng)建一個(gè)ID和數(shù)據(jù)結(jié)構(gòu)，將序號(hào)i作為基因的ID，每個(gè)基因設(shè)置一個(gè)狀態(tài)記錄該基因是否被訪(fǎng)問(wèn)，并且設(shè)置一個(gè)狀態(tài)記錄該基因是否為核心基因。每個(gè)基因的列表RGL記錄密度可達(dá)的鄰居基因，根據(jù)密度可達(dá)矩陣DRM計(jì)算RGL。設(shè)置一個(gè)狀態(tài)變量記錄該基因是否為邊緣基因，設(shè)置一個(gè)UID列表記錄子類(lèi)的UID值，將屬于子類(lèi)的基因添加到子類(lèi)對(duì)應(yīng)的UID列表中，如果基因的列表為空，說(shuō)明該基因尚未被分配子類(lèi)。

算法6是子聚類(lèi)算法的偽代碼。master節(jié)點(diǎn)將密度可達(dá)矩陣中的基因分配至各個(gè)worker節(jié)點(diǎn)，worker開(kāi)始并行迭代子聚類(lèi)程序。基因作為子聚類(lèi)程序的輸入，基因的子聚類(lèi)結(jié)果作為輸出。子聚類(lèi)程序?qū)⒑诵幕蜃鳛槌跏蓟宇?lèi)的中心基因Gi，對(duì)應(yīng)的UIDi將添加到核心基因Gi的UID列表中。如果一個(gè)密度可達(dá)基因Gf是核心基因，則將UIDi分配到Gf和Gf的所有密度可達(dá)基因。worker節(jié)點(diǎn)nextUID()函數(shù)生成唯一的UID，該函數(shù)保證不同worker程序UID值不相同。

算法6子聚類(lèi)算法

輸入：RDD。

//Spark RDD數(shù)據(jù)庫(kù)

輸出：UID。

//基因的子類(lèi)ID

1. for each (Gi?RDD&& CGstatus(Gi) == TRUE)

2. if (visiting_status(Gi) == UNVISITED)

//狀態(tài)未訪(fǎng)問(wèn)

3. visiting_status(Gi) == VISITED;

//狀態(tài)已訪(fǎng)問(wèn)

4. UID(i)←nextUID();

5. UIDlist(Gi) = UIDlist(Gi) + UID(i);

6. for eachGjin RGL_Gi

7. visiting_status(Gj) == VISITED;

8. UIDlist(Gk) = UIDlist(Gk) + UID(i);

9. if CGstatus(Gk) == TRUE

10.j=k；

11. end for

12. end if

13. end for

14. for each (Gi? RDD&&((UIDlist(Gi)== NULL)‖(visiting_status(Gj) == FALSE)))

15. BGstatus(Gi) = TRUE;

16. end for

每個(gè)子線(xiàn)程首先處理尚未訪(fǎng)問(wèn)的核心基因Gj，為Gj分配新的UIDj和密度可達(dá)基因集，UIDj的基因訪(fǎng)問(wèn)狀態(tài)變量修改為“訪(fǎng)問(wèn)”，最終訪(fǎng)問(wèn)完所有的基因。如果一個(gè)基因由兩個(gè)不同核心基因密度可達(dá)，則其UID列表中會(huì)存在多個(gè)UID值。

在程序結(jié)束階段將未被訪(fǎng)問(wèn)的基因標(biāo)記為邊緣基因BG。并且為每個(gè)子類(lèi)指定一個(gè)骨架基因，骨架基因的決定方法為同一個(gè)UID列表中所有基因的HFE值，計(jì)算方法為：

(14)

式中：Gpro為骨架基因;hGPro是子類(lèi)中所有基因的HFE;β是子類(lèi)的基因數(shù)量。

5.5 融合階段

將分布式處理獲得的子類(lèi)融合為總體的聚類(lèi)。計(jì)算兩個(gè)子類(lèi)UIDi、UIDj之間的相同成員，如果滿(mǎn)足以下任一條件：

(1) 子類(lèi)UIDi和UIDj之間相同基因的數(shù)量大于ξ，即{Ga∈UIDi,Ga∈UIDj,i≠j,|Ga|>ξ}；

(2) 子類(lèi)骨架基因之間的相似性大于τ，即{WρHFS(Gpro(UIDi),Gpro(UIDj))>τ,i≠j}。

那么將子類(lèi)UIDi和UIDj合并為新子類(lèi)UIDk。如果子類(lèi)UIDi和UIDj合并為UIDk，那么UIDi和UIDj的類(lèi)標(biāo)簽替換為UIDk，UIDk的基因數(shù)量等于UIDi和UIDj之和。

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

6.1 實(shí)驗(yàn)環(huán)境和參數(shù)設(shè)置

(1) 實(shí)驗(yàn)數(shù)據(jù)集和實(shí)驗(yàn)環(huán)境。在Hadoop YARN集群管理平臺(tái)搭建獨(dú)立的Apache Spark 2.0軟件。單一的基因微陣列數(shù)據(jù)集對(duì)實(shí)驗(yàn)的工具和環(huán)境具有敏感性，為了提高實(shí)驗(yàn)結(jié)果的可靠性，減小噪聲影響，采用多個(gè)微陣列數(shù)據(jù)集測(cè)試算法的性能。這些數(shù)據(jù)集組成了多條件、多實(shí)驗(yàn)下基因行為的猶豫模糊觀(guān)察樣本。

采用臨床數(shù)據(jù)和人工合成數(shù)據(jù)從多個(gè)角度測(cè)試算法的各個(gè)性能：臨床數(shù)據(jù)來(lái)自于不同的臨床實(shí)驗(yàn)，設(shè)立相似的實(shí)驗(yàn)條件以評(píng)價(jià)結(jié)果的精度和效果；人工合成數(shù)據(jù)為大規(guī)模數(shù)據(jù)集，加入噪聲數(shù)據(jù)測(cè)試算法的魯棒性和可擴(kuò)展性。第一個(gè)臨床數(shù)據(jù)集[23]為肝臟組織的cDNA、核苷酸、表達(dá)譜數(shù)據(jù)，數(shù)據(jù)集包含430個(gè)基因芯片微陣列、cDNA微陣列和寡核苷酸，將該數(shù)據(jù)集簡(jiǎn)稱(chēng)為肝臟。第二個(gè)臨床數(shù)據(jù)集[24]簡(jiǎn)稱(chēng)為POM，是體內(nèi)裂殖酵母全局細(xì)胞循環(huán)控制的基因微陣列數(shù)據(jù)集，該數(shù)據(jù)集是多組實(shí)驗(yàn)的微陣列數(shù)據(jù)，并且是一個(gè)關(guān)于時(shí)間的函數(shù)。第三個(gè)數(shù)據(jù)集[25]簡(jiǎn)稱(chēng)為SPC，是一個(gè)大規(guī)模的數(shù)據(jù)集，測(cè)試本算法的可擴(kuò)展性。

(2) RMA的參數(shù)設(shè)置。應(yīng)用輪盤(pán)賭策略從種群中選擇染色體，種群大小設(shè)為15，每對(duì)染色體應(yīng)用兩點(diǎn)交叉算子。染色體的第一個(gè)數(shù)量范圍為[n/2,3n/4]，n是種群的特征數(shù)量。采用精英機(jī)制保留父種群和子種群的最優(yōu)染色體，準(zhǔn)確率的權(quán)重設(shè)為特征數(shù)量權(quán)重的4倍，容錯(cuò)度α=2。在RMA優(yōu)化種群的過(guò)程中，保留種群中適應(yīng)度高于動(dòng)態(tài)閾值φ的i=3個(gè)最佳染色體。特征歸簡(jiǎn)處理之后，收集β%個(gè)最佳的特征組成新的特征空間。根據(jù)實(shí)驗(yàn)結(jié)果分析，將參數(shù)β設(shè)為定值5，最大迭代次數(shù)為20次，ρ=94.5，w1=4，w2=1，q=0.3，step=0.5，初始種群的特征數(shù)量范圍為[m/2, 3m/4]，m為染色體的特征數(shù)量。

(3) 聚類(lèi)處理的性能指標(biāo)。選擇3個(gè)常用的聚類(lèi)算法性能評(píng)價(jià)指標(biāo)，分別為： 鄧恩指數(shù)(Dunn Index,DI)、Davies-Bouldin指數(shù)(Davies-Bouldin Index,DBI)和Silhouette指數(shù)(Silhouette Index,SI)。

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

選擇4種較新的聚類(lèi)算法與本文算法做比較來(lái)綜合評(píng)價(jià)本文算法的性能，分別為MBC[25]、HRSC[26]、MOOC[27]、EGWOM[28]。其中MBC是基于建模的聚類(lèi)算法，HRSC是基于Hessian矩陣的聚類(lèi)算法，MOOC是基于多目標(biāo)優(yōu)化的聚類(lèi)算法，EGWOM則是基于灰狼優(yōu)化算法和MapReduce的聚類(lèi)算法。

(1) 聚類(lèi)性能結(jié)果。圖3、圖4、圖5分別為四個(gè)算法對(duì)于3個(gè)臨床數(shù)據(jù)集的實(shí)驗(yàn)結(jié)果。DBI值越小表示聚類(lèi)結(jié)果的簇內(nèi)緊密，簇間分離大。觀(guān)察實(shí)驗(yàn)結(jié)果可見(jiàn)，本文算法對(duì)于3個(gè)數(shù)據(jù)集的聚類(lèi)準(zhǔn)確率均優(yōu)于其他4種算法，本文算法在聚類(lèi)系統(tǒng)中采用猶豫模糊加權(quán)相關(guān)系數(shù)評(píng)估基因行為的相似性，根據(jù)基因的密度可達(dá)性將基因聚類(lèi)，該設(shè)計(jì)對(duì)于高維小樣本的基因微陣列數(shù)據(jù)實(shí)現(xiàn)了較好的效果。在特征歸簡(jiǎn)處理中，采用文化基因算法對(duì)特征空間進(jìn)行封裝式的特征選擇，實(shí)現(xiàn)了較好的優(yōu)化效果，降低了基因的冗余度且排除了噪聲基因。

圖3 肝臟數(shù)據(jù)集的聚類(lèi)性能結(jié)果

圖4 POM數(shù)據(jù)集的聚類(lèi)性能結(jié)果

圖5 MOOC數(shù)據(jù)集的聚類(lèi)性能結(jié)果

(2) 擴(kuò)展性實(shí)驗(yàn)和計(jì)算效率。采用POM數(shù)據(jù)集作為人工數(shù)據(jù)集，測(cè)試算法的可擴(kuò)展性和計(jì)算效率。將POM數(shù)據(jù)分別合成POM1、POM2、POM3、POM4四個(gè)數(shù)據(jù)集，POM1、POM2、POM3、POM4分別是將POM數(shù)據(jù)集復(fù)制1 000次、10 000次、100 000次和1 000 000次的數(shù)據(jù)集，結(jié)果如圖6所示?？梢钥闯?，隨著數(shù)據(jù)集規(guī)模10倍的擴(kuò)展，SI指標(biāo)并未呈現(xiàn)數(shù)倍的減低，數(shù)據(jù)規(guī)模增至1 000倍，SI指數(shù)大約降低0.1，DI指數(shù)大約降低0.15。結(jié)果證明了本文算法具有較好的可擴(kuò)展性，對(duì)于大規(guī)模的數(shù)據(jù)集依然具有較好的聚類(lèi)性能。

圖6 擴(kuò)展性實(shí)驗(yàn)的聚類(lèi)性能結(jié)果

圖7為不同規(guī)模數(shù)據(jù)集相應(yīng)的處理時(shí)間?？梢钥闯觯S著數(shù)據(jù)集規(guī)模的指數(shù)級(jí)增長(zhǎng)，處理時(shí)間也呈現(xiàn)出明顯的升高。比較POM1和POM4，數(shù)據(jù)集規(guī)模增加了1 000倍，但是時(shí)間僅提高了約10倍，所以本文算法在計(jì)算效率上也具有較好的擴(kuò)展性。

圖7 擴(kuò)展性實(shí)驗(yàn)的處理時(shí)間結(jié)果

(3) 與其他算法的處理效率比較。圖8為5個(gè)聚類(lèi)系統(tǒng)對(duì)于肝臟和SPC臨床數(shù)據(jù)集的處理時(shí)間結(jié)果?？梢钥闯觯琈BC和HRSC兩個(gè)算法的處理時(shí)間均明顯高于其他三個(gè)算法，主要原因是這兩個(gè)算法是一種串行處理機(jī)制，雖然MBC算法通過(guò)多模型化簡(jiǎn)減小了特征空間和數(shù)據(jù)維度，但是依然保持較高的計(jì)算時(shí)間。HRSC在Hessian矩陣的處理過(guò)程中需要大量的計(jì)算，具有較高的計(jì)算成本。MOOC是一種多目標(biāo)優(yōu)化的聚類(lèi)算法，該算法基于分布式多目標(biāo)框架加快了計(jì)算速度，EGWOM則基于MapReduce框架實(shí)現(xiàn)了分布式處理，但是MapReduce具有以下三點(diǎn)不足：① 每次迭代作為獨(dú)立任務(wù)重新處理，需要重新讀寫(xiě)數(shù)據(jù)和初始化數(shù)據(jù)；② 每次迭代需要重新從分布式文件系統(tǒng)加載程序，消耗I/O資源和CPU資源；③ 前一次迭代必須全部結(jié)束并且數(shù)據(jù)全部寫(xiě)入分布式文件系統(tǒng)，才能開(kāi)始下一次迭代,不僅具有較大的I/O開(kāi)銷(xiāo)，也導(dǎo)致了時(shí)間的延遲。本文算法基于Spark平臺(tái)實(shí)現(xiàn)了分布式的處理過(guò)程，對(duì)于迭代處理程序?qū)崿F(xiàn)了較快的速度。

圖8 5個(gè)聚類(lèi)系統(tǒng)的處理時(shí)間

7 結(jié) 語(yǔ)

為了支持海量高維小樣本數(shù)據(jù)的分析和挖掘，設(shè)計(jì)了基于遞歸文化基因和Spark分布式計(jì)算的高維大數(shù)據(jù)聚類(lèi)系統(tǒng)。該系統(tǒng)每次迭代作為獨(dú)立任務(wù)重新處理，不需要重新讀寫(xiě)數(shù)據(jù)和初始化數(shù)據(jù)；每次迭代直接從緩存系統(tǒng)加載程序，無(wú)須消耗額外的I/O資源和CPU資源。受益于遞歸文化基因算法和Spark分布式計(jì)算平臺(tái)的諸多優(yōu)點(diǎn)，本文算法實(shí)驗(yàn)獲得了較好的聚類(lèi)和聚類(lèi)效果。

Spark平臺(tái)將各個(gè)worker的數(shù)據(jù)緩存于內(nèi)存中，無(wú)須消耗額外的I/O資源和CPU資源，但是增加了內(nèi)存的負(fù)擔(dān)，在海量高維小樣本數(shù)據(jù)挖掘過(guò)程中產(chǎn)生大量的中間數(shù)據(jù)，導(dǎo)致不可忽略的內(nèi)存負(fù)擔(dān)，未來(lái)將關(guān)注于降低聚類(lèi)系統(tǒng)的內(nèi)存消耗。