CUDA-TP:基于GPU的自頂向下完整蛋白質(zhì)鑒定并行算法

段 瓊 田 博 陳 征 王 潔,2 何增有,2

1(大連理工大學(xué)軟件學(xué)院 遼寧大連 116620) 2 (遼寧省泛在網(wǎng)絡(luò)與服務(wù)軟件重點(diǎn)實(shí)驗(yàn)室(大連理工大學(xué)) 遼寧大連 116620) (wangjie1003@163.com)

Fig. 1 Comparison between TD and BU for protein identification圖1 蛋白質(zhì)鑒定中的TD與BU方法對(duì)比

蛋白質(zhì)組學(xué)(proteomics)是一門新興學(xué)科，主要研究細(xì)胞內(nèi)表達(dá)的所有蛋白質(zhì)及其活動(dòng)規(guī)律[1].由于基因的作用最終是由蛋白質(zhì)(protein)來(lái)體現(xiàn)，所以蛋白質(zhì)組學(xué)的研究有著更為重要的實(shí)際意義和價(jià)值.蛋白質(zhì)組學(xué)的一個(gè)重要目標(biāo)是能夠快速有效地進(jìn)行蛋白質(zhì)鑒定，即確定一個(gè)樣本中表達(dá)的所有的蛋白質(zhì).蛋白質(zhì)是由氨基酸分子鏈接而成的生物大分子，蛋白質(zhì)的一級(jí)結(jié)構(gòu)(氨基酸序列)唯一確定了蛋白質(zhì)的身份,因此可以通過(guò)識(shí)別氨基酸序列達(dá)到鑒定蛋白質(zhì)的目的.只有鑒定到生物樣品中真實(shí)表達(dá)的蛋白質(zhì)，才能準(zhǔn)確獲得蛋白質(zhì)相互作用、亞細(xì)胞定位等信息，為進(jìn)一步的疾病標(biāo)記物發(fā)現(xiàn)和新藥開(kāi)發(fā)提供強(qiáng)有力的支持[2].因此，蛋白質(zhì)鑒定是蛋白質(zhì)組學(xué)研究的基礎(chǔ)，對(duì)整個(gè)領(lǐng)域的進(jìn)一步發(fā)展和應(yīng)用有著十分重要的意義.

為解決蛋白質(zhì)鑒定問(wèn)題，目前以生物質(zhì)譜為基礎(chǔ)的蛋白質(zhì)組學(xué)分析方法主要有“自底向上”(bottom-up， BU)和“自頂向下”(top-down， TD)二種方法[3]，2種方法比對(duì)如圖1所示.BU方法通常先將蛋白質(zhì)的復(fù)雜樣品進(jìn)行酶切產(chǎn)生肽段的混合物；然后通過(guò)液相色譜等技術(shù)將這些肽段的混合物進(jìn)行分離，進(jìn)而通過(guò)質(zhì)譜技術(shù)將肽段碎裂，并根據(jù)碎裂譜圖中的離子峰信息進(jìn)行數(shù)據(jù)庫(kù)搜索來(lái)鑒定肽段；最后將鑒定到的肽段進(jìn)行組裝、推理獲得樣品中所含有的蛋白質(zhì)[4-6].

BU是由肽段推測(cè)蛋白質(zhì)序列，屬于“從局部推斷整體”.盡管BU已經(jīng)在當(dāng)前的蛋白質(zhì)組學(xué)研究中廣泛使用，但該類方法同樣存在著一系列的局限[7]：

1) 由于并不是所有的肽段都可以有效地被捕捉到并生成二級(jí)譜圖，導(dǎo)致很多蛋白質(zhì)沒(méi)有相應(yīng)的肽段鑒定結(jié)果，進(jìn)而無(wú)法鑒定到該蛋白質(zhì)的存在；

2) 即使可以通過(guò)少數(shù)幾個(gè)鑒定得到的肽段信息來(lái)推斷蛋白質(zhì)的存在與否，但是卻不能測(cè)繪出整個(gè)蛋白質(zhì)序列的全部信息，這主要包括蛋白質(zhì)各種氨基酸位點(diǎn)上的翻譯后修飾(post-translational modification, PTM)；

3) 由于是通過(guò)肽段的鑒定結(jié)果間接地得到蛋白質(zhì)的鑒定結(jié)果，而同一個(gè)肽段可以映射到多個(gè)不同的蛋白質(zhì)序列，這就需要解決一個(gè)復(fù)雜的計(jì)算問(wèn)題：蛋白質(zhì)推斷[8].而蛋白質(zhì)推斷問(wèn)題目前尚未徹底解決，仍有很多計(jì)算上的挑戰(zhàn)需要克服.

與BU策略相反，TD方法不需要將蛋白質(zhì)“切割”成更短的肽段序列，而是直接將完整的蛋白質(zhì)進(jìn)行分離和離子化，然后對(duì)其進(jìn)行碎裂并利用串聯(lián)質(zhì)譜測(cè)量由每個(gè)蛋白質(zhì)生成的二級(jí)譜圖.在數(shù)據(jù)分析階段，通過(guò)比對(duì)數(shù)據(jù)庫(kù)中的蛋白質(zhì)序列和實(shí)驗(yàn)二級(jí)譜圖進(jìn)行蛋白質(zhì)鑒定[9-10].這樣，通過(guò)完整蛋白質(zhì)的質(zhì)量及其碎裂譜的信息可以實(shí)現(xiàn)真正意義上的蛋白質(zhì)鑒定.另外，TD能夠保留多種PTM之間的關(guān)聯(lián)信息，逐漸成為蛋白質(zhì)組學(xué)研究中的熱點(diǎn)[5].

鑒于TD方法對(duì)于蛋白質(zhì)鑒定具有重要的生物學(xué)意義，針對(duì)TD策略中完整蛋白質(zhì)與譜圖的匹配問(wèn)題，已經(jīng)提出了一些有效的算法.在這些方法中，ProSight[11-12]采用泊松分布概率打分模型計(jì)算蛋白質(zhì)與譜圖的匹配概率，用“鳥(niǎo)槍標(biāo)注”的方法生成所有可能的蛋白質(zhì)變體數(shù)據(jù)庫(kù).這種方法最大的問(wèn)題是擴(kuò)展性差，如果變體增多可能導(dǎo)致“組合爆炸”.Frank等人[13]借鑒譜圖比對(duì)方法，以動(dòng)態(tài)規(guī)劃為基礎(chǔ)開(kāi)發(fā)了MS-TopDown算法，該算法能有效的鑒定蛋白質(zhì)及未知的PTM.Liu等人[14]在此基礎(chǔ)上通過(guò)優(yōu)化動(dòng)態(tài)規(guī)劃顯著提高了計(jì)算速度并提供了匹配的統(tǒng)計(jì)顯著性評(píng)估方法.在此后，TopPIC[15]，MASH suit[16]，MSpathFinder[17]，Ptop[18]等工具都是以動(dòng)態(tài)規(guī)劃為主體進(jìn)行蛋白質(zhì)數(shù)據(jù)庫(kù)搜索并在不同的方面做了改進(jìn).然而，其面臨的最大問(wèn)題還是匹配時(shí)間不盡如人意，這是由動(dòng)態(tài)規(guī)劃算法所固有的時(shí)間復(fù)雜度決定的.

高性能計(jì)算的發(fā)展對(duì)計(jì)算速度的提升做出了巨大貢獻(xiàn).自2003年開(kāi)始，圖形處理器(graphics proc-essing units, GPU)就在浮點(diǎn)運(yùn)算性能和存儲(chǔ)器帶寬上飛速發(fā)展.通用并行計(jì)算架構(gòu)(compute unified device architecture, CUDA)[19]是由英偉達(dá)公司推出的高性能運(yùn)算平臺(tái)，它能夠充分發(fā)揮GPU并行計(jì)算的優(yōu)勢(shì).隨著GPU可編程性的增強(qiáng)，GPU突破僅應(yīng)用于圖形處理領(lǐng)域的局限，開(kāi)始用于一些通用計(jì)算領(lǐng)域[20-21]，尤其是在生物信息學(xué)領(lǐng)域，為研究人員提供了新的思路.例如翟艷堂等人[22]提出的PTM鑒定算法MS-Alignment的并行架構(gòu)，有效地提高了在肽段上檢測(cè)PTM的效率；對(duì)于肽段鑒定問(wèn)題，Zhang等人[23]通過(guò)CUDA架構(gòu)將肽段鑒定算法RT-PSM[24]的速度提升了26倍；Baumgardner等人[25]成功對(duì)SpectraST[26]算法進(jìn)行了并行計(jì)算加速；Li等人[27]也順利地將譜圖點(diǎn)積算法并行化.近幾年來(lái)，CUDA并行計(jì)算架構(gòu)更是在基因序列比對(duì)方面取得了全新的進(jìn)展[28-30].

鑒于完整蛋白質(zhì)與質(zhì)譜數(shù)據(jù)匹配存在高耗時(shí)的缺點(diǎn)，本文借鑒譜圖比對(duì)的思想，以MS-TopDown算法為基礎(chǔ)，選擇CUDA編程模型設(shè)計(jì)CUDA-TP算法核心步驟，即完整蛋白質(zhì)與TD質(zhì)譜匹配分?jǐn)?shù)的計(jì)算.同時(shí)在計(jì)算步驟中引入平衡二叉搜索(adelson-velskii and landis， AVL)樹(shù)來(lái)保存每個(gè)路徑點(diǎn)上的信息并對(duì)MS-Filter[31]進(jìn)行改進(jìn)以加速蛋白質(zhì)過(guò)濾.最后，本文選取常用的target-decoy[32-33]方法對(duì)結(jié)果進(jìn)行錯(cuò)誤發(fā)現(xiàn)率(false discovery rate, FDR)控制.實(shí)驗(yàn)結(jié)果表明：在1%的FDR下CUDA-TP運(yùn)行速度分別是MS-TopDown和MS-Align+算法的10倍與2倍.

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

1) 提出了一種新型的基于GPU的完整蛋白質(zhì)與譜圖打分算法，該算法是對(duì)MS-TopDown的并行優(yōu)化，它繼承了可以鑒定未知PTM的特點(diǎn)，同時(shí)得到的最終分?jǐn)?shù)是譜圖網(wǎng)格中的全局分?jǐn)?shù)，并沒(méi)有使用MS-Align+方法中縮減譜圖網(wǎng)格搜索空間的策略.

2) 對(duì)MS-Filter算法進(jìn)行改進(jìn)加速了蛋白質(zhì)過(guò)濾過(guò)程.

3) 使用真實(shí)數(shù)據(jù)集進(jìn)行實(shí)驗(yàn)，驗(yàn)證了本文所提出的CUDA-TP算法的高效性.

1 問(wèn)題定義

1.1 符號(hào)表示

通常情況下，可以把一個(gè)譜圖轉(zhuǎn)換為其相應(yīng)的前綴殘基質(zhì)量(prefix residue mass， PRM)譜圖來(lái)表示完整蛋白質(zhì)的前綴肽段.例如二級(jí)串聯(lián)譜圖中一個(gè)質(zhì)量為ma的y離子單同位素譜峰在PRM譜圖中對(duì)應(yīng)的質(zhì)量為M-ma(M表示完整蛋白質(zhì)質(zhì)量).本文中，譜圖中所有譜峰都為單同位素譜峰，每個(gè)譜峰的單同位素質(zhì)量包含2個(gè)值：ma和其對(duì)應(yīng)的M-ma.同時(shí)，譜圖還包含1個(gè)空質(zhì)量0與代表完整蛋白質(zhì)PRM值的M-18(18表示1個(gè)水分子質(zhì)量).具體而言，可以把譜圖S表示為一個(gè)有序序列：

a0

其中,a0=0，an=M-18.相應(yīng)地，長(zhǎng)度為m的蛋白質(zhì)序列P表示為

b0

其中,bi表示蛋白質(zhì)序列P中從第1個(gè)氨基酸到第i個(gè)氨基酸的質(zhì)量之和(假設(shè)b0=0，bm=M′-18，M′為蛋白質(zhì)P未加任何修飾物的質(zhì)量).給定一個(gè)蛋白質(zhì)序列P和一個(gè)譜圖S，定義它們所組成的譜圖網(wǎng)格(spectra grid)為(a0,b0)，(a0,bm)，(an,b0)，(an,bm)四個(gè)點(diǎn)所組成的2維矩形方格，其共包含(n+1)(m+1)個(gè)點(diǎn).為了降低復(fù)雜性和方便計(jì)算，大部分的TD方法不包含譜峰強(qiáng)度，本文也將不考慮譜圖中譜峰強(qiáng)度.

1.2 質(zhì)譜匹配

蛋白質(zhì)序列P和譜圖S之間的匹配可以看作譜圖網(wǎng)格中的一條如圖2所示的路徑.路徑中一個(gè)點(diǎn)與其相鄰點(diǎn)的連接包含3種情況：

1) 如果從(ai,bj)到(ai′,bj′)滿足ai=ai′,bj=bj′，則呈斜對(duì)角線；

2) 如果從(ai,bj)到(ai′,bj′)滿足ai′>ai，則為垂直線，垂直線表示蛋白質(zhì)質(zhì)量增大，例如氨基酸位點(diǎn)上產(chǎn)生一個(gè)質(zhì)量大于0的PTM；

3) 如果從(ai,bj)到(ai′,bj′)滿足bi′>bi，則為水平線，水平線表示蛋白質(zhì)質(zhì)量減少，例如氨基酸位點(diǎn)上產(chǎn)生1個(gè)質(zhì)量小于0的PTM.匹配路徑上所經(jīng)過(guò)的點(diǎn)即為共享譜峰數(shù)量[34](shared peak count)，也稱為匹配分?jǐn)?shù).圖2(a)表示完美匹配，沒(méi)有發(fā)生任何PTM；圖2(b)表示第5和第6個(gè)氨基酸位點(diǎn)上分別產(chǎn)生-80Da和+90Da的PTM.

Fig. 2 Illustration of protein-spectra matching圖2 蛋白質(zhì)與譜圖之間的匹配

通過(guò)觀察可知，每條從(a0,b0)到(an,bm)的路徑都可能是蛋白質(zhì)序列P和譜圖S的一個(gè)有效匹配，而且匹配的路徑數(shù)目隨著質(zhì)量變化個(gè)數(shù)F(質(zhì)量變化對(duì)應(yīng)路徑中垂直與水平線)的增加而呈現(xiàn)指數(shù)增長(zhǎng).為了避免這種情況的發(fā)生，通常采用動(dòng)態(tài)規(guī)劃算法尋找譜圖網(wǎng)格中的最優(yōu)路徑[14-18]，這樣，算法的運(yùn)行時(shí)間只會(huì)隨著F的增加呈線性增長(zhǎng).

動(dòng)態(tài)規(guī)劃算法在尋找最優(yōu)路徑的過(guò)程中需要迭代地填充大小為n×m×F的數(shù)組D，D中的元素Di j(f)表示在最多允許發(fā)生f個(gè)質(zhì)量變化的前提下從(a0,b0)到(ai,bj)的最高匹配分?jǐn)?shù).數(shù)組D填充完成后，Dnm(F)的值即為蛋白質(zhì)序列P和譜圖S的匹配分?jǐn)?shù).在最優(yōu)路徑中，Di′j′(f′)的前驅(qū)點(diǎn)標(biāo)記為Di j(f)，如果2個(gè)數(shù)值對(duì)(ai,bj)和(ai′,bj′)滿足條件：

|ai-ai′-bj′+bj|<β，

(1)

則稱它們是余對(duì)角的(codiagonal)，其中β為誤差值.

通常，如果i

(2)

Di j(f)的計(jì)算為

Di j(f)=max(Ddiag(i,j)(f)+1，
Hi-1,j-1(f-1)+1).

(3)

Hi j(f)的狀態(tài)轉(zhuǎn)移方程為

Hi j(f)=max(Di j(f)，Hi-1,j(f),Hi,j-1(f)).

(4)

路徑起始點(diǎn)D0,0(f)=0.可以看出，質(zhì)量變化數(shù)F僅會(huì)讓數(shù)組D中的元素個(gè)數(shù)線性增加，算法的時(shí)間復(fù)雜度為O(nmF)，n和m分別表示譜圖S中譜峰數(shù)量與蛋白質(zhì)序列P的長(zhǎng)度.

為了減小匹配誤差，算法中的質(zhì)量變化可以選取已知的PTM(如雙氧化、磷酸化、甲基化)集合ΔPTM中的數(shù)值δ.因此，質(zhì)量變化分為了2部分：已知的PTM與未知的PTM.用Fs和Fg分別表示二者在匹配中所允許的最大個(gè)數(shù).至此，可以將diag(i,j)擴(kuò)展為diag(δ)(i,j)并使之滿足條件：

-β<|ai-ai′-bj′+bj|-δ<β，

(5)

則之前的數(shù)組D和H大小變?yōu)閚×m×Fg×Fs，Di j(g,s)中的值仍然表示從(0,0)到(i,j)的最大分?jǐn)?shù).把式(3)轉(zhuǎn)化成：

Di j(g,s)=max(Ddiag(i,j)(g,s)+1，
Hi-1,j-1(g-1,s)，Ddiag(δ)(i,j)(g,s-1)+1).

(6)

狀態(tài)轉(zhuǎn)移方程式(4)更新成

Hi j(g,s)=max(Di j(g,s)，
Hi-1,j(g,s),Hi,j-1(g,s)).

(7)

通過(guò)式(6)，可以利用動(dòng)態(tài)規(guī)劃算法，最終求得全局譜圖網(wǎng)格中完整蛋白質(zhì)序列P與譜圖S的匹配分?jǐn)?shù).

Fig. 3 The flowchart of CUDA-TP algorithm圖3 CUDA-TP算法流程圖

2 CUDA-TP算法

如第1節(jié)所述的算法雖然能在線性時(shí)間內(nèi)計(jì)算出匹配分?jǐn)?shù)，但當(dāng)?shù)鞍踪|(zhì)數(shù)據(jù)庫(kù)變大或者譜圖數(shù)量增多時(shí)，其時(shí)間開(kāi)銷依然不容樂(lè)觀.為了加速計(jì)算過(guò)程，本文提出了基于GPU的CUDA-TP算法，該算法首先在CPU端使用AVL樹(shù)優(yōu)化前驅(qū)節(jié)點(diǎn)的求取，然后設(shè)計(jì)并行模型加速式(6)的計(jì)算.同時(shí)，為了加快蛋白質(zhì)過(guò)濾，CUDA-TP還對(duì)蛋白質(zhì)過(guò)濾算法MS-Filter[31]做了改進(jìn).算法的整體流程如圖3所示.

2.1 計(jì)算diag

在譜圖網(wǎng)格中尋找最優(yōu)路徑時(shí)，首先要計(jì)算網(wǎng)格中每個(gè)點(diǎn)的前驅(qū)坐標(biāo)diag(i,j).由定義可知diag(i,j)必是(i,j)左上部余對(duì)角線上的點(diǎn).如圖4所示，為求A的前驅(qū)點(diǎn)，首先需要找到其左上部區(qū)域Z中與A在同一條余對(duì)角線上的點(diǎn)的集合{A0,A1,A2}，然后再在該集合內(nèi)選取diag(i,j)的坐標(biāo).在選取坐標(biāo)時(shí)，由式(6)可知，由于A1的diag值為A0的坐標(biāo)，所以匹配路徑中A1的匹配分?jǐn)?shù)(對(duì)應(yīng)D中元素)至少比A0大1，依次類推，得到點(diǎn)A的diag取值是離A最近點(diǎn)A2的坐標(biāo)值.因此，求譜圖網(wǎng)格中每個(gè)點(diǎn)的diag即為求距離該點(diǎn)最近且和它呈余對(duì)角的點(diǎn)的坐標(biāo).最簡(jiǎn)單的方法是遍歷滿足要求區(qū)域中的所有點(diǎn)，如圖4所示，A的遍歷區(qū)域?yàn)閆，但時(shí)間開(kāi)銷巨大.

Fig. 4 Calculate diag of all nodes in the spectra grid圖4 計(jì)算譜圖網(wǎng)格中所有點(diǎn)的diag

由式(1)易知，如果2個(gè)點(diǎn)是余對(duì)角的，那么它們所在坐標(biāo)對(duì)應(yīng)的數(shù)值對(duì)(ai,bj)的差值d在誤差范圍內(nèi)是相等的.也就是說(shuō)，集合{A0,A1,A2}對(duì)應(yīng)相同的d.利用這條性質(zhì)，可以通過(guò)AVL樹(shù)只遍歷一遍譜圖網(wǎng)格求取所有點(diǎn)的diag，如圖4所示.計(jì)算分步驟：

1) 將譜圖網(wǎng)格從左到右按列進(jìn)行劃分，每列元素的順序從上到下；

2) 建立一棵空的AVL樹(shù)，樹(shù)的節(jié)點(diǎn)包含1個(gè)坐標(biāo)和1個(gè)d值；

3) 從第1列的第1個(gè)元素開(kāi)始，依次計(jì)算坐標(biāo)元素的d值.如果AVL樹(shù)中沒(méi)有這個(gè)元素的d值，則把此值和它的坐標(biāo)作為新節(jié)點(diǎn)插入AVL樹(shù).如果查找到已經(jīng)有節(jié)點(diǎn)的d值與它相等，則把該元素的diag置為此節(jié)點(diǎn)所存儲(chǔ)的坐標(biāo)，同時(shí)更新節(jié)點(diǎn)的坐標(biāo)為該元素的坐標(biāo).

當(dāng)步驟3運(yùn)行至最后1列的最后1個(gè)元素時(shí)，所有點(diǎn)的diag求取完畢.由于AVL樹(shù)節(jié)點(diǎn)個(gè)數(shù)等于譜圖網(wǎng)格中所有節(jié)點(diǎn)的不同的d值數(shù)量，所以AVL樹(shù)的查找速度?O(lbnm)，能極大加快diag的計(jì)算.偽代碼如下：

算法1. 在CPU端計(jì)算diag.

輸入：譜圖a[n]、蛋白質(zhì)序列b[m]、誤差值β、AVL樹(shù)avl_tree；

輸出：所有點(diǎn)的diag.

①InitializeAVLtree(avl_tree)*將輸入的avl_tree初始化為null*

② fori=0 tomdo

③ forj=0 tondo

④d←b[i]-a[j]；

⑤ ifavl_treeincluded

⑥diag(i,j)←avl_tree[d](i,j)；

⑦avl_tree[d](i,j)←(i,j)；

⑧ else

⑨Insert(avl_tree,d,i,j)；

⑩ end if

需要注意的是，偽代碼中行⑤AVL樹(shù)在查找d值時(shí)，如果AVL樹(shù)節(jié)點(diǎn)中的值與該值的差的絕對(duì)值小于β，就認(rèn)為它們是相等的.avl_tree[d](i,j)表示AVL樹(shù)中查找到的等于d值的節(jié)點(diǎn)所包含的坐標(biāo).

2.2 計(jì)算diag(δ)

Fig. 5 Calculate diag(δ) of all nodes in the spectra grid圖5 計(jì)算譜圖網(wǎng)格中所有點(diǎn)的diag(δ)

不同于diag的計(jì)算，求解diag(δ)需要知道PTM集合ΔPTM中的所有數(shù)值δ.如果ΔPTM中有k個(gè)值，那么，求某個(gè)點(diǎn)的diag(δ)之前必須找出該點(diǎn)譜圖網(wǎng)格左上部與該點(diǎn)滿足式(5)的k個(gè)點(diǎn)集合，然后在這k個(gè)點(diǎn)集合中選取使匹配分值最大，即距離所求點(diǎn)最近的k個(gè)點(diǎn)的坐標(biāo)作為最終結(jié)果.如圖5所示，假設(shè)k=3，集合ΔPTM={δ1,δ2,δ3}，在求A點(diǎn)的diag(δ)時(shí)，首先在區(qū)域Y中尋找滿足條件的點(diǎn)集合{A0,A1,A2}，{B0,B1,B2}，{C0,C1,C2}，然后選取集合中距離A最近的3個(gè)點(diǎn)A2，B2，C2的坐標(biāo)作為返回值.顯然，算法1無(wú)法滿足這樣的計(jì)算要求.

由diag(δ)的計(jì)算過(guò)程可知，網(wǎng)格中的每個(gè)元素在求解時(shí)，其結(jié)果互不影響，即所有元素可以同時(shí)計(jì)算.這種性質(zhì)很好地符合了CUDA并行架構(gòu)要求.在CUDA架構(gòu)中，CPU作為主機(jī)(host)，GPU作為主機(jī)的協(xié)處理器(co-processor)，它被看作執(zhí)行高度線程化并行處理任務(wù)的計(jì)算設(shè)備(device).運(yùn)行在GPU上的CUDA并行計(jì)算函數(shù)稱為內(nèi)核(kernel)，以線程?hào)鸥?grid)形式組織，線程?hào)鸥裼扇舾蓚€(gè)線程塊(block)組成，每個(gè)線程塊又包含若干個(gè)線程(threads).實(shí)際運(yùn)行時(shí)，線程塊被分割成更小的線程束(wrap)來(lái)執(zhí)行運(yùn)算任務(wù)，1個(gè)線程束由連續(xù)的32個(gè)線程組成.

diag(δ)并行算法如圖5所示，網(wǎng)格中每個(gè)點(diǎn)的計(jì)算分配1個(gè)線程，線程之間并行執(zhí)行，線程塊與線程塊之間不需要同步，將所有的元素存入GPU的全局存儲(chǔ)器.對(duì)一個(gè)點(diǎn)左上部的搜索區(qū)域，仍然按照算法1進(jìn)行劃分.算法執(zhí)行過(guò)程中，始終維護(hù)1個(gè)包含k個(gè)點(diǎn)坐標(biāo)的數(shù)組g，該數(shù)組的元素為坐標(biāo)值且映射ΔPTM中的某個(gè)δ值.遍歷搜索區(qū)域時(shí)，如果2個(gè)點(diǎn)在同一個(gè)集合中，那么它們所對(duì)應(yīng)的δ必是相等的，則把后遍歷到的節(jié)點(diǎn)坐標(biāo)更新到數(shù)組g.如圖5中的點(diǎn)集合{A0,A1,A2}，它們都對(duì)應(yīng)δ2，當(dāng)遍歷至A1時(shí)，更新g[δ2]為A1的坐標(biāo)，依次類推，遍歷完成時(shí)，數(shù)組g即為A的最終結(jié)果.算法2描述了計(jì)算diag(δ)的具體過(guò)程.

算法2. 在GPU端計(jì)算diag(δ).

輸入：譜圖a[n]、蛋白質(zhì)序列b[m]、誤差值β、PTM的集合ΔPTM；

輸出：存儲(chǔ)所有點(diǎn)的diag(δ)數(shù)組g.

① (n_left,m_left)←GetCoord(blockDim.x,blockIdx.x,blockIdx.y)；*根據(jù)線程塊，線程信息獲取要處理的元素*

②d←fbs(b[m_left]-a[n_left])；

③ fori=0 tom_leftdo

④ forj=0 ton_leftdo

⑤delta←fbs(b[i]-a[j]-d)；

⑥ ifΔPTMincludedelta

⑦g[delta]←(i,j)；

⑧ end if

⑨ end for

⑩ end for

算法2中的行①表示從當(dāng)前的線程獲取所要計(jì)算的點(diǎn)坐標(biāo)；行②保存該坐標(biāo)對(duì)應(yīng)的d值；行③④開(kāi)始遍歷坐標(biāo)左上部區(qū)域的所有點(diǎn)；行⑤～求當(dāng)前遍歷點(diǎn)的delta值，如果PTM集合ΔPTM包含此值，則將其更新到數(shù)組g，最后返回的數(shù)組g即為最終結(jié)果；代碼的行⑥檢測(cè)ΔPTM中是否存在δ與delta之差的絕對(duì)值小于β，如果存在，則判斷集合ΔPTM包含delta.

2.3 數(shù)組計(jì)算并行架構(gòu)

為求蛋白質(zhì)序列P與譜圖S的匹配分?jǐn)?shù)，需要對(duì)式(6)和式(7)進(jìn)行計(jì)算，即求解數(shù)組D和數(shù)組H.在CUDA-TP算法中，D和H在邏輯上被聚集為數(shù)組E，E中每個(gè)元素包含2個(gè)值：Di,j，Hi,j.在計(jì)算Ei,j之前，需要得到Ei-1,j，Ei,j-1，Ei-1,j-1這3個(gè)值.如圖6所示，每次計(jì)算一個(gè)對(duì)角線上的元素，依次向下進(jìn)行，直至最后1個(gè)元素.不難發(fā)現(xiàn)，對(duì)角線上的元素計(jì)算時(shí)互相之間是并行的，某個(gè)對(duì)角線上的元素全部計(jì)算完成后，再進(jìn)行下一個(gè)對(duì)角線的計(jì)算.由此，可以設(shè)計(jì)并行架構(gòu)求解數(shù)組E.

Fig. 6 Execution of CUDA-TP algorithm圖6 CUDA-TP算法的執(zhí)行順序

若蛋白質(zhì)序列P的長(zhǎng)度為m，譜圖S的譜峰數(shù)量為n，則數(shù)組E的元素個(gè)數(shù)為nm.當(dāng)n或者m變得很大時(shí)，E的元素就會(huì)相應(yīng)的增多，導(dǎo)致要求解的空間很“龐大”，而GPU中同時(shí)并行執(zhí)行的線程數(shù)量是受限的，給每個(gè)元素指定1個(gè)線程顯然是不科學(xué)的.為了合理地利用GPU性能，可以把r×c個(gè)元素分為1個(gè)計(jì)算單元，每個(gè)單元分別分配1個(gè)線程，r和c根據(jù)GPU計(jì)算能力來(lái)指定大小(本文r和c均設(shè)置為2).線程按照標(biāo)號(hào)順序計(jì)算單元格中的元素，以保證在計(jì)算當(dāng)前元素時(shí)其左上部的所有元素是已知的.同時(shí)，并行的線程運(yùn)行在同一個(gè)線程塊內(nèi)，這樣能夠?qū)?shù)組E放到GPU的共享存儲(chǔ)器中.如圖7所示，計(jì)算單元A，B，C中的元素互不影響，它們是并行的，3個(gè)線程同時(shí)對(duì)其進(jìn)行計(jì)算且每個(gè)線程按順序計(jì)算4個(gè)元素.

Fig. 7 Parallel processing architecture圖7 并行計(jì)算架構(gòu)

算法3. 在GPU端計(jì)算數(shù)組E.

輸入：譜圖a[n]、蛋白質(zhì)序列b[m]、計(jì)算單元大小r和c，diag，diag(δ)；

輸出：所求數(shù)組E.

①idx←GetIdx(blockDim.x,blockIdx.x,blockIdx.y)；

②p←min(nr,mc)；

③ fori=0 tordo

④ forj=0 tocdo

⑤ (ei,ej)←GetEIndex(E,p,idx,i,j)；

⑥ if (ei,ej)≠?

⑦E[ei,ej]←GetMax(E,ei,ej,diag,diag(δ))；*獲取結(jié)果并填充進(jìn)數(shù)組E*

⑧ end if

⑨ end for

⑩ end for

算法3的行①首先獲取當(dāng)前的線程標(biāo)號(hào)；行②求得最大并行數(shù)p；行③④表示線程按順序計(jì)算r×c個(gè)元素；行⑤得到要處理的元素坐標(biāo)，如果該坐標(biāo)代表的元素不為空元素，則將數(shù)組E中此元素的數(shù)值更新；最后的行返回?cái)?shù)組E.需要注意的是，更新數(shù)組元素值的函數(shù)GetMax中的參數(shù)diag為算法1的結(jié)果，diag_delta即為3.2節(jié)所述的diag(δ)，E包含的數(shù)組D和H值嚴(yán)格按照式(6)和式(7)求解.

2.4 時(shí)間復(fù)雜度分析

數(shù)組填充時(shí)，串行程序依次計(jì)算數(shù)組E中的元素，完成F個(gè)E數(shù)組的計(jì)算需要O(NF).由2.3節(jié)可知，CUDA-TP算法將E劃分為F(r×c)個(gè)單元格，每次有p個(gè)線程同時(shí)運(yùn)行且每個(gè)線程在單元格中串行執(zhí)行，但由于空元素的存在，平均同時(shí)運(yùn)行線程的數(shù)量約為p2, GPU執(zhí)行完畢后最多剩余(r-1)n+(c-1)m個(gè)元素，于是CUDA-TP計(jì)算E的時(shí)間復(fù)雜度為

(8)

因此，質(zhì)譜匹配的總時(shí)間復(fù)雜度由O(N2+FN)變?yōu)?/p>

(9)

2.5 蛋白質(zhì)過(guò)濾

對(duì)譜圖S來(lái)說(shuō)，要從蛋白質(zhì)數(shù)據(jù)庫(kù)中尋找1個(gè)蛋白質(zhì)P，使得它與這個(gè)蛋白質(zhì)的匹配分?jǐn)?shù)最高.蛋白質(zhì)數(shù)據(jù)庫(kù)通常包含很多個(gè)蛋白質(zhì)，讓所有蛋白質(zhì)與譜圖S進(jìn)行匹配分值計(jì)算，時(shí)間開(kāi)銷是巨大的.常用的方法是為每個(gè)譜圖挑選出L個(gè)候選蛋白質(zhì)，然后譜圖與這L個(gè)候選蛋白質(zhì)計(jì)算匹配分?jǐn)?shù)并以最高分值的蛋白質(zhì)作為鑒定結(jié)果.MS-Filter[31]是目前進(jìn)行蛋白質(zhì)過(guò)濾很有效的方法(蛋白質(zhì)過(guò)濾問(wèn)題定義及相關(guān)算法參見(jiàn)文獻(xiàn)[31])，它為每個(gè)蛋白質(zhì)與譜圖S計(jì)算過(guò)濾分?jǐn)?shù)，以分?jǐn)?shù)高低作為是否過(guò)濾蛋白質(zhì)的條件.CUDA-TP通過(guò)優(yōu)化MS-Filter算法3個(gè)計(jì)算步驟中的步驟2來(lái)加速蛋白質(zhì)過(guò)濾，優(yōu)化后的計(jì)算過(guò)程：

1) 在低精度下計(jì)算譜圖與蛋白質(zhì)的卷積數(shù)組；

2) 查找對(duì)角線分值大于閾值的候選區(qū)域，如果相鄰候選區(qū)域的最大邊界與最小邊界值之差小于所有氨基酸中甘氨酸的最小質(zhì)量75.05，則將這2個(gè)區(qū)域合并；

3) 在高精度下依次對(duì)合并后的候選區(qū)域計(jì)算卷積數(shù)組，如果卷積數(shù)組中對(duì)角線分值大于設(shè)定的閾值，則將閾值更新為此分值，最后的閾值即為蛋白質(zhì)的過(guò)濾分?jǐn)?shù).

蛋白質(zhì)的過(guò)濾加速體現(xiàn)在步驟2的候選區(qū)域合并，通過(guò)合并分散的小區(qū)域來(lái)減少步驟3高精度卷積數(shù)組的個(gè)數(shù).

CUDA-TP創(chuàng)建L個(gè)流(stream)對(duì)過(guò)濾出的候選蛋白質(zhì)與譜圖的匹配分?jǐn)?shù)并行計(jì)算，CUDA架構(gòu)的流并行是粗粒度的并行，它能夠使GPU運(yùn)算的同時(shí)與CPU進(jìn)行數(shù)據(jù)交換.如圖8所示，為每個(gè)候選蛋白質(zhì)分配1個(gè)流，流與流之間并行執(zhí)行.

Fig. 8 Stream parallelism for calculating the matching score between the spectrum and candidate proteins圖8 譜圖與候選蛋白質(zhì)匹配分?jǐn)?shù)計(jì)算的流并行

3 實(shí)驗(yàn)與結(jié)果

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

本文采用人類細(xì)胞蛋白p65[35]、大腸桿菌蛋白[36](escherichia coli， Ecoli)、人類核心組蛋白(human core histones)H4[37]和鼠傷寒沙門氏菌[38](salmonella typhimurium， ST)質(zhì)譜數(shù)據(jù)來(lái)進(jìn)行實(shí)驗(yàn)，所有的蛋白質(zhì)數(shù)據(jù)庫(kù)均來(lái)自美國(guó)國(guó)立生物信息中心NCBI①.原始質(zhì)譜數(shù)據(jù)文件使用ReAdw②和MS-Deconv[39]轉(zhuǎn)化為只包含單同位素譜峰的譜圖.

由于TD方法是對(duì)完整蛋白質(zhì)的鑒定，譜峰數(shù)量太少會(huì)導(dǎo)致匹配分?jǐn)?shù)過(guò)低，所以只保留前體質(zhì)量(precursor mass)大于2500Da且至少包含10個(gè)譜峰的譜圖[14].最終的實(shí)驗(yàn)數(shù)據(jù)集及其分布如表1和圖9所示.

Table1TheNumberofSpectrafromFourDatasetsandtheNumberofProteinsintheCorrespondingDatabases

表1 數(shù)據(jù)集中的譜圖數(shù)量及數(shù)據(jù)庫(kù)中蛋白質(zhì)個(gè)數(shù)

Fig. 9 Distribution of the number of peaks and the length of proteins圖9 譜峰數(shù)量和蛋白質(zhì)長(zhǎng)度分布

p65數(shù)據(jù)集包含240個(gè)譜圖，Ecoli數(shù)據(jù)集包含921個(gè)譜圖，H4和ST分別包含1 245，4 339個(gè)譜圖，譜圖中的譜峰數(shù)量分布如圖9(a)所示.p65與Ecoli數(shù)據(jù)集的譜峰較為均勻地分布在20～190之間;ST數(shù)據(jù)集的譜峰數(shù)量大多在20～80之間，占比57.2%，為2 482個(gè)；而H4數(shù)據(jù)集的譜峰數(shù)量多數(shù)分布在160～200之間，占比57.8%，為720個(gè).

4個(gè)實(shí)驗(yàn)數(shù)據(jù)集的蛋白質(zhì)數(shù)據(jù)庫(kù)分別含有349，2 206，2 433，4 606個(gè)蛋白質(zhì)，蛋白質(zhì)序列長(zhǎng)度分布如圖9(b)所示.p65蛋白質(zhì)數(shù)據(jù)庫(kù)中的蛋白質(zhì)長(zhǎng)度大多分布在240～480之間，占比60.2%，共211個(gè)；ST蛋白質(zhì)數(shù)據(jù)庫(kù)中的蛋白質(zhì)長(zhǎng)度在40～240之間的有2 947個(gè)，占整個(gè)數(shù)據(jù)庫(kù)的64.1%；H4和Ecoli蛋白質(zhì)數(shù)據(jù)庫(kù)中的蛋白質(zhì)長(zhǎng)度多集中在80～340，H4包含1 989個(gè)，Ecoli包含1 883個(gè)，各自占比為81.4%，83.3%.

3.2 CUDA-TP運(yùn)行時(shí)間

CUDA-TP基于譜圖比對(duì)思想，與MS -Align+③通過(guò)減少搜索譜圖網(wǎng)格空間以達(dá)到運(yùn)行時(shí)間降低的策略不同，它是串行算法MS-TopDown的并行加速算法.本文通過(guò)設(shè)置不同的參數(shù)F(常用策略是固定Fg，改變Fs)和L來(lái)對(duì)CUDA-TP與MS-TopDown的運(yùn)行時(shí)間進(jìn)行比對(duì)，實(shí)驗(yàn)環(huán)境如表2所示:

Table 2 The Running Environment表2 實(shí)驗(yàn)環(huán)境

Fig. 10 Running time of CUDA-TP and MS-TopDown on four datasets with different parameters 圖10 4個(gè)數(shù)據(jù)集上CUDA-TP與MS-TopDown在不同參數(shù)下的運(yùn)行時(shí)間

MS-TopDown沒(méi)有蛋白質(zhì)過(guò)濾步驟，選用2.5節(jié)改進(jìn)的MS-Filter算法進(jìn)行蛋白質(zhì)過(guò)濾.MS-Align+從官方網(wǎng)站下載，MS-Align+只提供參數(shù)F的設(shè)置，運(yùn)行時(shí)L為軟件中的默認(rèn)參數(shù)，因此只對(duì)比不同F(xiàn)參數(shù)下的運(yùn)行時(shí)間.同時(shí)，實(shí)驗(yàn)采用常用的target-decoy[32-33]進(jìn)行FDR控制，所得實(shí)驗(yàn)結(jié)果FDR值均小于1%.本文首先設(shè)定不同的允許最多發(fā)生的PTM個(gè)數(shù)F和蛋白質(zhì)過(guò)濾數(shù)量L，從多角度對(duì)比CUDA-TP與MS-TopDown的運(yùn)行時(shí)間，4個(gè)數(shù)據(jù)集的實(shí)驗(yàn)結(jié)果如圖10所示:

在p65數(shù)據(jù)集上，CUDA-TP平均運(yùn)行時(shí)間為22 min，MS-TopDown平均運(yùn)行時(shí)間為212 min，CUDA-TP速度是MS-TopDown的9.6倍；在Ecoli數(shù)據(jù)集上，CUDA-TP平均運(yùn)行時(shí)間為110 min，MS-TopDown平均運(yùn)行時(shí)間為1 212 min，CUDA-TP速度是MS-TopDown的11倍；在ST數(shù)據(jù)集上，CUDA-TP平均運(yùn)行時(shí)間為243 min，MS-TopDown平均運(yùn)行時(shí)間為2 315 min，CUDA-TP速度是MS-TopDown的9.5倍；在H4數(shù)據(jù)集上，CUDA-TP平均用時(shí)139 min，MS-TopDown用時(shí)1 405 min，CUDA-TP速度是MS-TopDown的10.1倍.可以看出，雖然質(zhì)譜數(shù)據(jù)集ST是H4的將近4倍，但是用時(shí)卻只是H4的2倍，這是由于H4譜圖中譜峰數(shù)量與其數(shù)據(jù)庫(kù)蛋白質(zhì)長(zhǎng)度普遍比ST的高，導(dǎo)致要求解的譜圖網(wǎng)格元素增多，增加了運(yùn)行時(shí)間.

MS-Align+只可以調(diào)整F的大小，L默認(rèn)固定為20，因此本文對(duì)比分析了3種方法在不同F(xiàn)參數(shù)下的運(yùn)行時(shí)間，實(shí)驗(yàn)結(jié)果如圖11所示:

Fig. 11 Running time of three methods on four datasets with different parameters圖11 4個(gè)數(shù)據(jù)集上3種方法在不同參數(shù)下的運(yùn)行時(shí)間

3種方法的平均運(yùn)行時(shí)間在表3中給出，從表3中可以看到CUDA-TP的運(yùn)行速度約是MS-Align+的2倍，這證明基于譜圖比對(duì)思想的并行計(jì)算方法明顯優(yōu)于通過(guò)減少譜圖網(wǎng)格搜索空間來(lái)加速計(jì)算過(guò)程的策略.

Table 3 Average Running Time of Three Methods表3 3種方法平均運(yùn)行時(shí)間 min

3.3 CUDA-TP算法吞吐率

并行算法除了運(yùn)行時(shí)間，有時(shí)更加關(guān)心其運(yùn)行時(shí)的吞吐率.MS-TopDown與MS-Align+是單核CPU 程序，并沒(méi)有進(jìn)行多線程優(yōu)化，這是由于2種方法采用以空間換取時(shí)間的策略，單核運(yùn)行時(shí)就幾乎達(dá)到了電腦資源的占用極限，實(shí)驗(yàn)中二者平均內(nèi)存使用率在90%以上，多核CPU的運(yùn)行幾乎是不可完成的.文獻(xiàn)[18]也指出單核MS-Align+運(yùn)行大規(guī)模人類蛋白質(zhì)質(zhì)譜數(shù)據(jù)時(shí)，內(nèi)存需求甚至高達(dá)40 GB.而CUDA-TP并行算法可以在顯存1 GB的電腦上流暢運(yùn)行，其內(nèi)存占用不超過(guò)4 GB，因此具有更廣的適用性.

CUDA-TP的時(shí)間復(fù)雜度已經(jīng)在2.4節(jié)詳細(xì)給出，3種方法在不同數(shù)據(jù)集上的吞吐率如表4所示，吞吐率指單位時(shí)間內(nèi)算法執(zhí)行的計(jì)算單元數(shù)量.從表4中可以看出CUDA-TP吞吐率約為MS-TopDown和MS-Align+的11倍，這在Ecoli和H4數(shù)據(jù)上表現(xiàn)尤為明顯，說(shuō)明譜圖網(wǎng)格元素的增多雖然增加了運(yùn)行時(shí)間，但吞吐率并沒(méi)有隨著降低.MS-Align+通過(guò)減少求取網(wǎng)格的數(shù)量來(lái)減少計(jì)算時(shí)間，但吞吐率并沒(méi)有實(shí)質(zhì)提升.因此，CUDA-TP在運(yùn)行時(shí)間和算法的吞吐率上要明顯優(yōu)于目前的MS-TopDown和MS-Align+方法.

Table 4 Throughput of Three Methods 表4 3種方法的吞吐率 10-6cells

Table 4 Throughput of Three Methods 表4 3種方法的吞吐率 10-6cells

DatasetCUDA-TPMS-TopDownMS-Align+p65108.59.89.7Ecoli120.410.511.3H4115.111.211.7ST106.78.99.2

MS-TopDown在數(shù)據(jù)集上的執(zhí)行時(shí)間通常要花費(fèi)大量時(shí)間(本文中MS-TopDown在最大數(shù)據(jù)集上的運(yùn)行時(shí)間在2 d左右)，這大大降低了其在蛋白質(zhì)鑒定中的實(shí)用性，而本文提出的基于GPU的蛋白質(zhì)鑒定并行算法運(yùn)行速度是其10倍，有效地解決了譜圖比對(duì)思想應(yīng)用于大規(guī)模數(shù)據(jù)的弊端，與現(xiàn)有的MS-Align+算法相比具有明顯優(yōu)勢(shì)，通過(guò)上述多種不同的實(shí)驗(yàn)測(cè)試，驗(yàn)證了CUDA-TP的優(yōu)秀性能.CUDA-TP源代碼托管在github公共服務(wù)網(wǎng)站①.

4 總結(jié)及展望

當(dāng)前，“自頂向下”的蛋白質(zhì)組學(xué)飛速發(fā)展，已經(jīng)成為大規(guī)模鑒定蛋白質(zhì)及定位PTM的關(guān)鍵技術(shù)，但這些技術(shù)的應(yīng)用算法在運(yùn)行時(shí)間上還存在瓶頸.本文研究了TD策略下的蛋白質(zhì)鑒定問(wèn)題，提出了一種新型的基于GPU的完整蛋白質(zhì)鑒定并行算法CUDA-TP.1)該算法通過(guò)流并行和優(yōu)化MS-Filter來(lái)加速蛋白質(zhì)過(guò)濾；2)引入AVL樹(shù)降低譜圖網(wǎng)格中每個(gè)元素前驅(qū)節(jié)點(diǎn)的求取時(shí)間；3)在GPU端設(shè)計(jì)了計(jì)算迭代公式的CUDA并行架構(gòu).實(shí)驗(yàn)結(jié)果表明CUDA-TP可以有效地加速完整蛋白質(zhì)鑒定，與通過(guò)減少譜圖搜索空間來(lái)?yè)Q取效率的MS-Align+相比具有明顯優(yōu)勢(shì).

在TD策略下對(duì)完整蛋白質(zhì)進(jìn)行鑒定，除了計(jì)算蛋白質(zhì)與譜圖的匹配分?jǐn)?shù)，還需要進(jìn)一步評(píng)估匹配結(jié)果的統(tǒng)計(jì)顯著性.因此如何計(jì)算匹配分?jǐn)?shù)的同時(shí)得到蛋白質(zhì)與譜圖匹配的統(tǒng)計(jì)顯著性評(píng)估結(jié)果，并且不降低時(shí)間運(yùn)行效率，這將是我們下一步的研究工作.