生物效應大數(shù)據(jù)評估聚類算法的并行優(yōu)化

彭紹亮，楊順云，孫哲，程敏霞，崔英博，王曉偉，李非，伯曉晨，廖湘科

1. 湖南大學信息科學與工程學院&國家超級計算長沙中心, 湖南 長沙 410082；

2. 國防科技大學計算機學院, 湖南 長沙 410073；

3. 中國人民解放軍軍事醫(yī)學科學院, 北京 100850

1 引言

近年來，隨著生物技術的發(fā)展，生物信息的數(shù)據(jù)量達到了一個更高的級別，生物醫(yī)藥領域的實驗手段和研究方法均發(fā)生了巨大的變革，呈現(xiàn)出“大數(shù)據(jù)”的趨勢，傳統(tǒng)的單機計算已經(jīng)不足以應對海量的數(shù)據(jù)和繁重的計算任務。對于大數(shù)據(jù)處理，常用的思路是并行計算，其包括多進程和多線程兩種并行等級。生物效應分析流程主要包括比對和聚類。本文主要針對大量藥物化合物制劑刺激下人體細胞反應的基因表達譜數(shù)據(jù)，完成細胞反應大數(shù)據(jù)的分析處理。主要分為以下3個步驟。

● 數(shù)據(jù)預處理：利用開源工具1KTools對整合網(wǎng)絡細胞印記庫（library of integrated network-based cellular signatures，LINCS）的原始基因譜數(shù)據(jù)進行預處理，得到實驗核心程序能夠使用的數(shù)據(jù)格式并寫出文件。

● 基因探針富集分析（gene set enrichment analysis，GSEA）算法的核心實現(xiàn)：利用預處理后的數(shù)據(jù)完成富集積分矩陣的計算，采用MPI+OpenMP二級并行的策略負載均衡地劃分數(shù)據(jù)，充分利用資源完成計算，并按進程寫出結果文件。

● 并行聚類：以比對結果為輸入，實現(xiàn)K-medoids[1]聚類算法及其優(yōu)化，并對每次迭代過程同樣利用MPI+OpenMP二級并行的策略進行并行化加速，最后將聚類結果寫出到文件，每個表達譜歸屬于某一聚類。

2 相關工作

2.1 生物效應評估方法

隨著生物技術的飛速發(fā)展，特別是以新一代測序技術為代表的高通量分析技術的發(fā)展，生命科學的年數(shù)據(jù)產(chǎn)出能力已經(jīng)達到PB級，呈現(xiàn)出“大數(shù)據(jù)”的趨勢，涉及海量的組學數(shù)據(jù)、文獻數(shù)據(jù)、臨床數(shù)據(jù)等。僅公開的數(shù)據(jù)庫（如GEO[2]、ArrayExpress[3]、TCGA[4]等）就包含了大量病原微生物感染刺激下人體細胞反應的基因表達譜數(shù)據(jù)。2010年美國國立衛(wèi)生研究院（NIH）啟動了LINCS項目[5]，其目標是系統(tǒng)地檢測15 000種化學分子對15種典型人體細胞刺激后的基因表達情況。目前該計劃第一期已獲得15種典型細胞中3 000余個基因沉默和5 000余種化學小分子刺激下的130余萬個全基因組表達譜。

“生物效應評估”字面上理解就是評估這些生物制劑對人體細胞產(chǎn)生的效應，具體而言就是指評估這種生物制劑究竟會致使某種疾病還是治愈某種疾病。從而僅通過計算手段快速確定相關的檢測標識物和治療靶標，極大地縮短防治手段的研發(fā)過程，以快速有效地應對可能的生物威脅，給人類健康提供更多的保障。

對于轉(zhuǎn)錄組數(shù)據(jù)的比較指標，采用了GSEA[6-8]算法中提出的富集積分，它基于排序的Kolmogorov-Smirnov檢驗統(tǒng)計量計算方法，并且采用顯著性分析、多重假設檢驗的方法對得到的富集積分進行統(tǒng)計分析，衡量結果的可靠性。

2.2 高性能計算技術在生物效應評估中的需求

目前GSEA在表達譜分析中得到廣泛應用，隨著 RNA-seq和低成本轉(zhuǎn)錄組L1000技術的流行，越來越多的大規(guī)模轉(zhuǎn)錄組數(shù)據(jù)出現(xiàn)，對于這樣大規(guī)模的數(shù)據(jù)分析研究，往往需要快速的GSEA計算過程以支持數(shù)據(jù)挖掘和機器學習應用。于是，為了應對大數(shù)據(jù)場景下快速計算的需求，就有了利用超級計算對計算過程進行分布式并行加速的必要。

目前國內(nèi)高性能計算技術也取得了豐碩的成果，其與世界先進水平高性能計算機之間的差距正在逐步縮小，6次蟬聯(lián)超級計算機Top 500第一名[9]的“天河二號”代表著我國超級計算機的卓越成績。

3 算法介紹

3.1 GSEA算法

GSEA算法主要用于分析兩個不同表形樣本集之間的表達差異，其基本思想是檢驗所定義基因集（gene set）S中的基因在整個微陣列實驗測得的已排序的所有基因列表（gene list）L中是均勻分布的還是集中于頂端或底部的。其本身是非常豐富、全面的，包含大量的統(tǒng)計學分析手段，分為3個主要步驟：富集積分的計算、估計效應量（effect size，ES）顯著性水平、調(diào)整多重假設檢驗。

其核心是使用Kolmogorov-Smirnov統(tǒng)計量進行富集積分計算，反映某個基因集在整個排序列表的頂部或底部集中出現(xiàn)的程度。ES的計算是在基因列表L中順序步移的，從初始值ES(S)開始，當步移中遇到S中的基因時，增加ES(S)，否則，減小ES(S)。增加或減小的幅度依賴于基因與表型間的相關程度。ES就是ES(S)整個步移過程中與0的最大偏差的絕對值最大的值。

ES的計算過程[8]如圖1所示。

富集積分的計算框架總結如下。

根據(jù)樣本數(shù)據(jù)表達相關性，對含有N個基因的表達譜進行排序，得到有序的基因譜L={g1,g2,…,gn}，使得序列第j個位置的基因表達量就是第j個基因gj的表達量，即有式r(gj)=rj成立。

計算基因集（基因探針）S在序列L上的hit向量和miss向量，計算式如式（1）和式（2）所示。

當隨機選取基因集S時，ES(S)會相對較小，但是如果其聚集在L的尾部或者頂端，會有一個很大的ES(S)值。當p=0（指數(shù)）的時候，ES(S)將會減變?yōu)闃藴蔏olmogorov-Smirnov統(tǒng)計分布（這是一個判別兩個未知分布的總體是否有同一分布的非參數(shù)檢驗）。

在實際實現(xiàn)時，預排序工作和計算Phit、Pmiss按照式（1）和式（2）計算。而找最大偏離處的ES值，一般是先對Phit和Pmiss兩個向量進行前綴求和，求和的過程中再對當前迭代次的前綴項作差，記錄最大的差值項。當遍歷完整個序列L時，就得到了最后的富集積分ES。

基于求解富集積分的算法描述過程可以清晰地分析算法各部分的時間復雜度，假設序列L長度為n，S的長度為m，算法的完成主要包含以下3個步驟的工作。

步驟1 對原始基因譜進行表達量排序：一般是達到O（nlogn）。

步驟2 計算Phit、Pmiss向量：每判斷一個表達譜基因是否命中就要掃描整個探針S，故而時間復雜度為O（mn）。

步驟3 求解ES：算法描述中也提到需要對步驟2中兩個向量進行遍歷計算前綴和，故而時間復雜度為O（n）。

綜上所述，GSEA是對表達譜進行分析的有效手段，但目前已有的實現(xiàn)工具（如R語言實現(xiàn)的GSEA-P-R.1.0[10]，Python實現(xiàn)的SAM-GS[11]以及Matlab實現(xiàn)的GSEA2[12]）都受限于腳本語言的低效性，難以達到需求的計算速度，同時它們對原本算法的實現(xiàn)極其直接，沒有進行任何的性能優(yōu)化，也沒有使用任何的并行化加速手段，這就導致其對海量表達譜數(shù)據(jù)的分析難以滿足實際的時間需求。因而急需更加有效的并行加速手段完成GSEA的快速分析[13]。

3.2 K-medoids聚類算法

經(jīng)典的K-means[14]算法處理大數(shù)據(jù)集合時非常高效且伸縮性較好。但當聚類的樣本中有“噪聲”（離群點）時，會產(chǎn)生較大的誤差，“噪聲”對確定新一輪質(zhì)心影響較大，造成所得質(zhì)點和實際質(zhì)點位置偏差過大。為了解決該問題，本文采用的K-medoids[1]聚類算法提出了新的質(zhì)點選取方式，而不是像K-means算法一樣簡單地采用均值計算法。

在K-medoids算法中，每次迭代后的質(zhì)點都是從聚類的樣本點中選取的，而選取的標準就是當該樣本點成為新的質(zhì)點后能提高類簇的聚類質(zhì)量，使類簇更緊湊[15]。該算法使用絕對誤差標準來定義一個類簇的緊湊程度。而在實現(xiàn)中是直接選取當前類簇中與其他元素平均距離最近的點作為新的聚類中心。另一方面，之所以選擇這個聚類算法也是因為它始終從已有的點中尋找新的中心，這就意味著計算過程中不會產(chǎn)生新的中心點，也就無需重新計算相似度，進而也就不需要返回比對算法部分重新計算富集積分。

然而，該算法具有同K-means算法一樣的一些缺陷，比如：K-medoids也需要隨機地產(chǎn)生初始聚類中心，不同的初始聚類中心可能導致完全不同的聚類結果。本文通過K-medoids++算法來優(yōu)化這部分計算。

圖1 ES的計算過程

K-medoids++算法選擇初始質(zhì)心的基本思想是：初始的聚類中心之間的相互距離要盡可能地遠。其算法描述如下：

步驟1 從輸入的數(shù)據(jù)點集合中隨機選擇一個點作為第一個聚類中心；

步驟2 對于數(shù)據(jù)集中的每一個點x，計算它與最近聚類中心（指已選擇的聚類中心）的距離D(x)；

步驟3 選擇一個新的數(shù)據(jù)點作為新的聚類中心，D(x)較大的點，被選取作為聚類中心的概率較大；

步驟4 重復步驟2和步驟3直到k個聚類中心被選出來；

步驟5 利用這k個初始聚類中心來運行標準K-medoids算法。

為了分析其復雜性，首先必須明確由于算法的執(zhí)行存在隨機性，不能形式化地判斷它經(jīng)過多少次迭代后收斂，因此對于算法復雜度的分析只針對單次迭代過程，假設有n個數(shù)據(jù)點和k個中心，復雜度分析如下。

● 生成初始聚類中心：對于每個數(shù)據(jù)點都要遍歷已有的每一個聚類中心，從而找到離它最近的中心，然后從這幾個最近中心中找到那個距離最遠的點作為新的初始中心，如此重復k次。初略估計復雜度約為：O（k（kn+n））=O（k2n）（如果不是K-medoids++，直接隨機生成這一步基本沒有開銷）。

● 劃分數(shù)據(jù)到類簇：每個點遍歷各聚類中心O（kn）。

● 尋找新的聚類中心：每個點要在各自的類簇中計算它到其他點的平均距離，然后再在各自的類簇中確定平均距離最小的點作為新的聚類中心，綜合來看這一步的時間復雜度是O（n2）。

綜上所述，開銷最大的是尋找新中心的部分，至于生成初始中心的過程，因為始終只有一次，當?shù)螖?shù)比較多的時候可以忽略不計。

4 基于MPI和OpenMP的并行優(yōu)化

4.1 基因譜比對算法的優(yōu)化

針對標準富集積分的計算在并行和循環(huán)過程中可能存在的優(yōu)化部分，本文在成功移植GSEA算法后，有針對性地進行了優(yōu)化，主要分為優(yōu)化單獨的計算例程和消除多次計算中出現(xiàn)的冗余計算兩部分。

4.1.1 優(yōu)化富集積分標準計算例程

由于實驗過程會重復地計算富集積分，優(yōu)化該步驟勢必會獲得較為理想的性能加速效果，所以不再根據(jù)GSEA算法按部就班地直接實現(xiàn)富集積分的計算，而是采用下面的優(yōu)化策略。

（1）只考察命中位置

回顧富集積分的計算過程，通過在基因列表L中順序步移，從初始值ES(S)開始，當步移中遇到S中的基因時，則增加ES(S)，否則，減小ES(S)。增加或減小的幅度依賴于基因與表型間的相關程度。ES就是ES(S)整個步移過程中與0的最大偏差的絕對值最大的值。通過觀察分析不難發(fā)現(xiàn)，最大偏差位置只會出現(xiàn)在命中位置附近。若命中基因的前一個基因未被命中，則其前一個基因有可能為低峰；若其后一個基因未被命中，則其后一位置可能為目前的高峰。

（2）對命中基因排序

命中向量的命中位置對于計算富集積分至關重要，在計算開始前，對命中向量根據(jù)其命中位置排序可以直接得出命中位置的ES(S)值，即其未被命中的個數(shù)為其位置減1，進一步省去標準計算例程中的前綴求和部分。通過優(yōu)化，可以把GSEA中富集積分前綴求和的復雜度由O（n）降低為O（lbm）+O（m）。而在實際情況中，基因譜的長度一般來說遠大于上下調(diào)基因集的長度，也就是n遠大于m。當大規(guī)模重復進行富集積分標準計算流程時，這樣的優(yōu)化將會帶來十分可觀的性能提升。

4.1.2 消除冗余計算

（1）預排序

因為任務是基因譜的兩兩比對，所以同一個基因譜在計算過程中肯定會被重復用到。而對同一個基因譜的排序工作也會因此反復地進行，這些都是沒有必要的工作。于是在讀入文件后就先對所有的基因譜進行預排序，之后處理的都將是排序后的轉(zhuǎn)錄組基因譜，從而排除冗余的排序過程。

（2）建立位置索引

富集積分的計算過程首先計算一個命中向量，即一個基因譜的上調(diào)或下調(diào)基因集在另一個排好序的基因譜中出現(xiàn)的位置，而這步操作需要循環(huán)遍歷基因譜和基因集。如果基因集的長度是m，基因譜的長度是n，這樣該步驟將是一個時間復雜度為O（mn）的操作。為了提高效率，系統(tǒng)先掃描一遍基因譜建立每個基因的位置索引數(shù)組，然后再掃描一遍基因集即可完成工作，從而將時間復雜度減小到O（m+n），并且用索引數(shù)組替代原來的排序基因譜也并不會造成額外的空間開銷。

（3）構建三元組保存基因譜結構

與排序一樣的問題，同一個基因譜會因為兩兩比對反復地建立索引，這還是冗余的操作。同樣地，在完成預排序的同時就先為每個基因譜建立索引并用之替代原來的排序基因譜。但是只有索引數(shù)組并不能確定原來基因的上下調(diào)基因集，故而在讀入數(shù)據(jù)后的預處理部分，用一個三元組保存基因譜的結構，它分別由上調(diào)基因集、下調(diào)基因集和索引數(shù)組三部分組成。

4.2 基因譜比對算法的并行優(yōu)化

本實驗的并行并不是對單次計算富集積分的算法過程進行的，而是通過有效的數(shù)據(jù)劃分和負載均衡手段對大規(guī)模計算富集積分矩陣的過程進行的。其原因是：第一，單次富集積分的算法本身不適合并行，雖然有些前綴求和的操作也能通過消除循環(huán)依賴的辦法強行并行，但是這會增加整體的工作量，得不償失；第二，只計算一個富集積分即使在全基因長度2萬多的基因譜上也不是個很大的工作量，對它并行粒度太小，反而會大量增加額外調(diào)度開銷。

為了達到計算過程負載均衡的目標，實驗首先對數(shù)據(jù)在進程間進行合理的劃分。因為軟件的輸入是兩個基因譜集的文件，所以數(shù)據(jù)的劃分其實就是對這兩個文件的劃分，使每個進程擁有大致等量的待計算基因譜。

基于此，對于文件1，系統(tǒng)直接按進程數(shù)進行劃分，使每個進程持有文件1的基因譜數(shù)相差不超過1，如果文件1的基因譜數(shù)剛好能夠被啟動的進程數(shù)整除，則它將被均勻地劃分給每個進程，這是容易做到的。對于文件2，如果再將它按進程進行負載均衡的劃分，將不能完成兩文件中任意兩個基因譜的比對工作。比如分給進程0的文件1的數(shù)據(jù)將不能與分給進程1的文件2的數(shù)據(jù)進行比對，如果強行比對，各進程在計算的過程中還要進行大規(guī)模的通信工作，這樣做的開銷是巨大的。于是，在一般內(nèi)存足夠的情況下，選擇犧牲一定空間的策略，讓每個進程持有文件2的全部基因譜數(shù)據(jù)，最大限度地保障了系統(tǒng)的計算性能。

相應地，在進程內(nèi)部，采用多線程的策略對文件2的數(shù)據(jù)進行負載均衡的劃分與計算。因為線程任務先天共享內(nèi)存，這樣，就可以充分發(fā)揮多核并行的優(yōu)勢，完成大規(guī)模并行計算工作。

整個數(shù)據(jù)劃分方式如圖2所示。

圖2清楚地顯示出兩個文件的基因譜數(shù)據(jù)在各進程以及進程內(nèi)部的線程中的劃分方式，同時可以看到，最后每個進程會將結果寫到各個子矩陣行。

對于具體的實現(xiàn)，文件1直接用封裝好的I/O函數(shù)定位相應的部分數(shù)據(jù)進行讀取。對于文件2，提供了3種通信方式。第一種方式是直接由每個進程讀取文件2的全部數(shù)據(jù)，這樣將不存在任何通信工作，但是也沒有體現(xiàn)任何的并行工作。第二種方式是讓一個進程讀取文件2的全部數(shù)據(jù)，然后將之均勻地劃分給每個進程，但這樣不僅要像前一種方法一樣等待一個進程讀完一個完整文件2，還要再進行通信，顯然數(shù)據(jù)劃分的性能并不理想。第三種方式是，讓每個進程一開始只并行地讀取文件2的部分數(shù)據(jù)，然后通過全局通信操作讓每個進程持有文件2的全部數(shù)據(jù)，這樣讀文件的時間將被大大地縮短。如果全局操作的實現(xiàn)足夠高效，即使加上額外的通信，也將獲得可觀的性能提升。在實驗中可以首先對比3種通信方式的運行效率，然后選擇最佳方式繼續(xù)實驗數(shù)據(jù)的測試工作。

4.3 聚類算法的并行實現(xiàn)

圖2 總體數(shù)據(jù)劃分

聚類實驗的并行化就是對K-medoids算法的并行化，沒有太多的優(yōu)化技巧，只在比對結果的基礎上先由每個進程讀入自己的那部分ES矩陣行，每一行代表一個基因譜相對于其他所有基因譜的距離向量，在后面的算法過程中，每個進程都只進行自己這幾行基因譜的計算，其中會用集合通信的方式在各進程間全局維護一個類標記向量，這就是利用信息傳遞接口（MPI）完成的進程級并行。至于每個進程對持有的基因譜集進行劃分以充分利用單節(jié)點處理器資源完成類簇規(guī)劃和尋找新中心的操作，是OpenMP完成的線程級并行工作。每次聚類迭代并行化實現(xiàn)框架如圖3所示。

如圖3所示，在并行實現(xiàn)的過程中仍有許多細節(jié)需要注意，其具體實現(xiàn)如下所述。

步驟1 每個進程讀取其下部分ES矩陣結果（由其序號，故而這部分進程數(shù)應該與之前計算ES矩陣并行寫文件時一致）。

圖3 每次聚類迭代并行化實現(xiàn)框架

步驟2 每個進程下生成一個長度為n1的類標記向量local_classflag，作為全局類標記向量的局部，所有進程的向量總長應為基因譜總數(shù)n。同時每個進程內(nèi)應該保有自己每行基因譜的global_rank起始號，由基因譜總數(shù)與進程數(shù)就可計算判定。

步驟3 生成0～n的k個不重復的隨機數(shù)，作為k個初始聚類中心。

步驟4 劃分類：每個進程利用OpenMP并行加速，判斷進程中每一行表達譜的歸屬類（根據(jù)其到k個聚類中心的距離，選取距離最小的聚類中心作為歸屬類），將local_classflag對應位置標記為該聚類中心編號。

步驟5 找新聚類中心：合并local_classflag為global_classflag到0號進程，然后廣播給其他進程（當然如果是平均劃分的可以直接Allgather）。每個進程將local_avedis對應的位置設置成其每一行基因譜到同一歸屬類其他基因譜的平均距離，然后將這些局部平均長度向量合并到0號進程的global_avedis向量，其長度為基因譜總數(shù)n。找到每類中平均長度最小的基因譜作為新的聚類中心。

步驟6 判斷和之前相比，聚類中心是否發(fā)生變化，若無變化則停止，并輸出相應結果；否則，重復執(zhí)行步驟4。

值得注意的是，這里基因譜之間的距離是用富集積分的倒數(shù)來進行衡量的，ES值越大說明基因譜越相似，距離就越近，故而這樣處理很容易理解。對于優(yōu)化的K-Medoids++算法，只是在步驟3中進行更多的操作，按算法介紹的流程生成初始聚類中心。

5 實驗結果分析

實驗的原始表達譜數(shù)據(jù)來自LINCS項目，現(xiàn)已提供在GEO官網(wǎng)。

它有在各種實驗條件下得到的多種規(guī)模的表達譜數(shù)據(jù)，一般是以.gct或者.gctx為后綴的HDF5文件格式。因為實驗中直接采用1KTools開源工具進行原始數(shù)據(jù)解析，所以并不關注該文件本身的復雜結構，只需解析所需的表達譜數(shù)據(jù)即可。lKTools項目現(xiàn)在還在維護之中，提供了R、Java、Matlab和Python 4種語言的解析包，本文主要使用的是Matlab版本。

程序中主要解析了表達譜的3類數(shù)據(jù)。

● mat：基因譜集的表型矩陣。

● rid：每個基因的標識。

● cid：每個表達譜的標識。

5.1 基因譜兩兩比對，計算富集積分矩陣實驗分析

實驗中在“天河二號”上使用了包含2萬基因譜的數(shù)據(jù)集進行實驗，設置不同的并行程度，記錄各部分的運行時間，以研究程序的性能和可拓展性。實驗結果見表1。

由表1可知，數(shù)據(jù)載入部分的時間開銷是比較小的，說明并行文件讀入還是比較高效的。另外，對3種通信模式的性能進行分析，可以看到，全局通信優(yōu)于無通信且優(yōu)于點對點通信。分析3種模式實現(xiàn)策略不難發(fā)現(xiàn)這樣的結果還是比較合理的。首先點對點通信效果最差，是由于它要先等0號進程讀完整個文件后，再由0號進程將數(shù)據(jù)發(fā)給其他進程。因為在無通信的情況下，它只花費讀文件的時間，所以這種點對點通信策略的性能會次于無通信策略，但是它能減少系統(tǒng)總體的通信開銷。

全局通信每個進程只讀取部分文件，然后通過Allgather操作將數(shù)據(jù)整合并發(fā)給全部進程，這樣整體的I/O是比較小的，并且可以并行完成，只是多了大量的通信。顯然它是有可能在性能上優(yōu)于無通信策略的。而測試的結果也確實如此。

表1 基因譜比對實驗各階段運行記錄表

寫文件的操作以進程數(shù)為基準劃分，可以看到，在相同進程數(shù)下，寫文件的時間是大致相同的。而不同進程數(shù)下，測試結果也表現(xiàn)出比較好的可拓展性。

ES矩陣的計算部分也體現(xiàn)出較好的可拓展性。前期看來，基本是每增加一倍的并行度，運行時間就縮短一半，效率接近于1，理想上已經(jīng)趨近于強可拓展的應用。整體來看，效果比較理想。當然，不能保證繼續(xù)增加并行度，這樣的效果還可以繼續(xù)保持。所以擴大數(shù)據(jù)規(guī)模和并行程度，用5萬的基因譜集進行對比分析，結果如圖4所示。

從圖4中可以看出，數(shù)據(jù)集規(guī)模越大，在擴大并行度時，并行效率保持得越高。這意味著在集群規(guī)模足夠大的情況下，本文實驗的程序可極好地適應大規(guī)模數(shù)據(jù)分析的需求。

5.2 基因譜聚類實驗分析

（1）性能分析

使用2萬規(guī)模和5萬規(guī)模基因譜數(shù)據(jù)集比較單次迭代過程中兩種聚類算法的執(zhí)行效率，如圖5所示。

因為算法的隨機性，相同參數(shù)下多次運行程序執(zhí)行的收斂步數(shù)是不一樣的，所以不能用總的執(zhí)行時間評估算法的性能，只能如圖5一樣根據(jù)單次迭代執(zhí)行時間來分析，并將其轉(zhuǎn)換成并行效率。

不難看出，在每次迭代的過程中算法都保持了比較良好的可拓展性，前期算法都會維持接近于1的高效率。同時在數(shù)據(jù)集規(guī)模相同的情況下，K-medoids和K-medoids++的實現(xiàn)效率也基本接近，它們只在尋找初始聚類中心時不同，并不會對后面每次的迭代過程產(chǎn)生太大影響。

（2）算法收斂性評估

有效的聚類算法必須能夠快速地收斂，該部分實驗在60核的條件下對不同聚類算法的收斂步數(shù)與執(zhí)行時間進行評估，其結果如圖6所示?？梢园l(fā)現(xiàn)聚類數(shù)越多，收斂越慢，執(zhí)行時間也越長。另外，由于K-medoids++在生成初始聚類中心時盡可能地保證了樣本空間距離足夠遠，所以它相對于K-medoids能夠更快地收斂，效果也更佳。同時，也能夠發(fā)現(xiàn)在相同聚類數(shù)目下，數(shù)據(jù)規(guī)模越大，收斂也可能越快。

另外值得注意的是，如果算法在非相鄰的兩次迭代中得到了相同的新聚類中心，將導致其不收斂的情況發(fā)生。為了避免出現(xiàn)這種情況，改進的實驗中維護了一個聚類中心集列表，每次迭代中產(chǎn)生的未出現(xiàn)過的新聚類中心將被添加其中，如果新中心已經(jīng)存在于列表中，則算法收斂。但該策略的缺點也是十分明顯的：首先，迭代多次后，該列表會越來越長，消耗大量內(nèi)存；同時，遍歷列表判斷算法是否收斂的開銷也會越來越大。因此建議即使數(shù)據(jù)集足夠大，也不要把聚類中心數(shù)設置得太多（比如超過500個），這樣算法還是能夠在開銷過大之前有效地收斂的。

圖4 不同規(guī)模數(shù)據(jù)集并行效率比較

圖5 單次迭代聚類算法并行效率比較

圖6 不同聚類中心數(shù)目下收斂步數(shù)和執(zhí)行時間

6 結束語

本文通過比對和聚類分析實現(xiàn)生物效應的快速評估，通過優(yōu)化和并行加速，在極短的時間內(nèi)完成細胞反應大數(shù)據(jù)的分析處理。同時，對比2萬數(shù)據(jù)量和5萬數(shù)據(jù)量的結果可以發(fā)現(xiàn)，在比對和聚類算法的實現(xiàn)中，擴大并行度時，數(shù)據(jù)量越大，并行效率越高且保持得越久，算法收斂越快，展現(xiàn)了實驗的可拓展性和收斂性?？傮w來說，實驗達到了預期目標。

[1]PARK H S, JUN C H. A simple and fast algorithm for K-medoids clustering[J].Expert Systems with Applications, 2009,36(2): 3336-3341.

[2]BARRETT T, TROUP D B, WILHITE S E,et al. NCBI GEO: archive for functional genomics data sets[J]. Nucleic Acids Research,2013, 41(Database Issue): 991-995.

[3]PARKINSON H, KAPUSHESKY M,SHOJATALAB M, et al. ArrayExpress-a public database of microarray experiments and gene expression profiles[J]. Nucleic Acids Research, 2007, 35(Database Issue): 747-750.

[4]KATARZYNA T, PATRYCJA C, MACIEJ W. The Cancer Genome Atlas (TCGA):an immeasurable source of knowledge[J].Contemporary Oncology, 2015, 19(1A):68-77.

[5]WON S J, WU H C, LIN K T, et al.Discovery of molecular mechanisms of lignan justicidin a using L1000 gene expression profiles and the Library of Integrated Network-based Cellular Signatures database[J]. Journal of Functional Foods, 2015, 16: 81-93.

[6]張威, 張揚, 曹文君, 等. GAGE和GSEA在基因集研究中的有效性比較[J]. 現(xiàn)代生物醫(yī)學進展, 2013, 13(10): 1849-1852.ZHANG W, ZHANG Y, CAO W J, et al.Comparative study of GAGE and GSEA in Gene-set analysis[J]. Progress in Modern Biomedicine, 2013, 13(10): 1849-1852.

[7]馮春瓊, 鄒亞光, 周其趙, 等. GSEA在全基因組表達譜芯片數(shù)據(jù)分析中的應用[J]. 現(xiàn)代生物醫(yī)學進展, 2009, 9(13): 2553-2557.FENG C Q, ZOU Y G, ZHOU Q Z, et al.The application of GSEA in data analysis of Genome microarray[J]. Progress in Modern Biomedicine, 2009, 9(13): 2553-2557.

[8]SUBRAMANIAN A, TAMAYO P,MOOTHA V K, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles[J]. Proceedings of the National Academy of Sciences of the United States of America, 2005, 102(43): 15545-15550.

[9]DONGARRA J. Visit to the National University for Defense Technology Changsha, China[R]. 2013.

[10]MOOTHA V K, LINDGREN C M,ERIKSSON K F, et al. PGC-1alpharesponsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes[J].Nature Genetics, 2003, 34(3): 267-273.

[11]DINU I, POTTER J D, MUELLER T, et al.Improving gene set analysis of microarray data by SAM-GS[J]. BMC Bioinformatics,2007, 8(1): 1-13.

[12]CARRO M S, WEI K L, ALVAREZ M J,et al. The transcriptional network for mesenchymal transformation of brain tumors[J]. Nature, 2010, 463(7279): 318-325.

[13]GAGGERO M, LEO S, MANCA S, et al.Parallelizing bioinformatics applications with MapReduce[J]. ResearchGate, 2008:1-6.

[14]MACQUEEN J B. Some Methods for classification and Analysis of Multivariate Observations[C]//The 5th Berkeley Symposium on Mathematical Statistics and Probability, June 21-July 18, 1965,California, USA. Berkeley: University of California Press, 1967: 281-297.

[15]劉小鳳. 適用于RNA二級結構預測的改進K-medoids聚類算法研究[D]. 秦皇島: 燕山大學, 2014.LIU X F. Research of suitableK-medoids clustering algorithm for RNA secondary structures prediction[D]. Qinhuangdao:Yanshan University, 2014.