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

1. 湖南大學(xué)信息科學(xué)與工程學(xué)院&國(guó)家超級(jí)計(jì)算長(zhǎng)沙中心, 湖南 長(zhǎng)沙 410082；

2. 國(guó)防科技大學(xué)計(jì)算機(jī)學(xué)院, 湖南 長(zhǎng)沙 410073；

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

1 引言

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

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

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

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

2 相關(guān)工作

2.1 生物效應(yīng)評(píng)估方法

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

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

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

2.2 高性能計(jì)算技術(shù)在生物效應(yīng)評(píng)估中的需求

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

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

3 算法介紹

3.1 GSEA算法

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

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

ES的計(jì)算過(guò)程[8]如圖1所示。

富集積分的計(jì)算框架總結(jié)如下。

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

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

當(dāng)隨機(jī)選取基因集S時(shí)，ES(S)會(huì)相對(duì)較小，但是如果其聚集在L的尾部或者頂端，會(huì)有一個(gè)很大的ES(S)值。當(dāng)p=0（指數(shù)）的時(shí)候，ES(S)將會(huì)減變?yōu)闃?biāo)準(zhǔn)Kolmogorov-Smirnov統(tǒng)計(jì)分布（這是一個(gè)判別兩個(gè)未知分布的總體是否有同一分布的非參數(shù)檢驗(yàn)）。

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

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

步驟1 對(duì)原始基因譜進(jìn)行表達(dá)量排序：一般是達(dá)到O（nlogn）。

步驟2 計(jì)算Phit、Pmiss向量：每判斷一個(gè)表達(dá)譜基因是否命中就要掃描整個(gè)探針S，故而時(shí)間復(fù)雜度為O（mn）。

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

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

3.2 K-medoids聚類算法

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

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

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

圖1 ES的計(jì)算過(guò)程

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

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

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

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

步驟4 重復(fù)步驟2和步驟3直到k個(gè)聚類中心被選出來(lái)；

步驟5 利用這k個(gè)初始聚類中心來(lái)運(yùn)行標(biāo)準(zhǔn)K-medoids算法。

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

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

● 劃分?jǐn)?shù)據(jù)到類簇：每個(gè)點(diǎn)遍歷各聚類中心O（kn）。

● 尋找新的聚類中心：每個(gè)點(diǎn)要在各自的類簇中計(jì)算它到其他點(diǎn)的平均距離，然后再在各自的類簇中確定平均距離最小的點(diǎn)作為新的聚類中心，綜合來(lái)看這一步的時(shí)間復(fù)雜度是O（n2）。

綜上所述，開(kāi)銷(xiāo)最大的是尋找新中心的部分，至于生成初始中心的過(guò)程，因?yàn)槭冀K只有一次，當(dāng)?shù)螖?shù)比較多的時(shí)候可以忽略不計(jì)。

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

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

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

4.1.1 優(yōu)化富集積分標(biāo)準(zhǔn)計(jì)算例程

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

（1）只考察命中位置

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

（2）對(duì)命中基因排序

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

4.1.2 消除冗余計(jì)算

（1）預(yù)排序

因?yàn)槿蝿?wù)是基因譜的兩兩比對(duì)，所以同一個(gè)基因譜在計(jì)算過(guò)程中肯定會(huì)被重復(fù)用到。而對(duì)同一個(gè)基因譜的排序工作也會(huì)因此反復(fù)地進(jìn)行，這些都是沒(méi)有必要的工作。于是在讀入文件后就先對(duì)所有的基因譜進(jìn)行預(yù)排序，之后處理的都將是排序后的轉(zhuǎn)錄組基因譜，從而排除冗余的排序過(guò)程。

（2）建立位置索引

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

（3）構(gòu)建三元組保存基因譜結(jié)構(gòu)

與排序一樣的問(wèn)題，同一個(gè)基因譜會(huì)因?yàn)閮蓛杀葘?duì)反復(fù)地建立索引，這還是冗余的操作。同樣地，在完成預(yù)排序的同時(shí)就先為每個(gè)基因譜建立索引并用之替代原來(lái)的排序基因譜。但是只有索引數(shù)組并不能確定原來(lái)基因的上下調(diào)基因集，故而在讀入數(shù)據(jù)后的預(yù)處理部分，用一個(gè)三元組保存基因譜的結(jié)構(gòu)，它分別由上調(diào)基因集、下調(diào)基因集和索引數(shù)組三部分組成。

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

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

為了達(dá)到計(jì)算過(guò)程負(fù)載均衡的目標(biāo)，實(shí)驗(yàn)首先對(duì)數(shù)據(jù)在進(jìn)程間進(jìn)行合理的劃分。因?yàn)檐浖妮斎胧莾蓚€(gè)基因譜集的文件，所以數(shù)據(jù)的劃分其實(shí)就是對(duì)這兩個(gè)文件的劃分，使每個(gè)進(jìn)程擁有大致等量的待計(jì)算基因譜。

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

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

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

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

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

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

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

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

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

步驟1 每個(gè)進(jìn)程讀取其下部分ES矩陣結(jié)果（由其序號(hào)，故而這部分進(jìn)程數(shù)應(yīng)該與之前計(jì)算ES矩陣并行寫(xiě)文件時(shí)一致）。

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

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

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

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

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

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

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

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

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

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

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

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

● rid：每個(gè)基因的標(biāo)識(shí)。

● cid：每個(gè)表達(dá)譜的標(biāo)識(shí)。

5.1 基因譜兩兩比對(duì)，計(jì)算富集積分矩陣實(shí)驗(yàn)分析

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

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

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

表1 基因譜比對(duì)實(shí)驗(yàn)各階段運(yùn)行記錄表

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

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

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

5.2 基因譜聚類實(shí)驗(yàn)分析

（1）性能分析

使用2萬(wàn)規(guī)模和5萬(wàn)規(guī)?；蜃V數(shù)據(jù)集比較單次迭代過(guò)程中兩種聚類算法的執(zhí)行效率，如圖5所示。

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

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

（2）算法收斂性評(píng)估

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

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

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

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

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

6 結(jié)束語(yǔ)

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

[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]張威, 張揚(yáng), 曹文君, 等. GAGE和GSEA在基因集研究中的有效性比較[J]. 現(xiàn)代生物醫(yī)學(xué)進(jìn)展, 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在全基因組表達(dá)譜芯片數(shù)據(jù)分析中的應(yīng)用[J]. 現(xiàn)代生物醫(yī)學(xué)進(jìn)展, 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二級(jí)結(jié)構(gòu)預(yù)測(cè)的改進(jìn)K-medoids聚類算法研究[D]. 秦皇島: 燕山大學(xué), 2014.LIU X F. Research of suitableK-medoids clustering algorithm for RNA secondary structures prediction[D]. Qinhuangdao:Yanshan University, 2014.