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

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

1(大連理工大學軟件學院 遼寧大連 116620) 2 (遼寧省泛在網絡與服務軟件重點實驗室(大連理工大學) 遼寧大連 116620) (wangjie1003@163.com)

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

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

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

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

1) 由于并不是所有的肽段都可以有效地被捕捉到并生成二級譜圖，導致很多蛋白質沒有相應的肽段鑒定結果，進而無法鑒定到該蛋白質的存在；

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

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

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

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

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

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

本文的主要貢獻有3個方面：

1) 提出了一種新型的基于GPU的完整蛋白質與譜圖打分算法，該算法是對MS-TopDown的并行優化，它繼承了可以鑒定未知PTM的特點，同時得到的最終分數是譜圖網格中的全局分數，并沒有使用MS-Align+方法中縮減譜圖網格搜索空間的策略.

2) 對MS-Filter算法進行改進加速了蛋白質過濾過程.

3) 使用真實數據集進行實驗，驗證了本文所提出的CUDA-TP算法的高效性.

1 問題定義

1.1 符號表示

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

a0

其中,a0=0，an=M-18.相應地，長度為m的蛋白質序列P表示為

b0

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

1.2 質譜匹配

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

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

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

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

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

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

動態規劃算法在尋找最優路徑的過程中需要迭代地填充大小為n×m×F的數組D，D中的元素Di j(f)表示在最多允許發生f個質量變化的前提下從(a0,b0)到(ai,bj)的最高匹配分數.數組D填充完成后，Dnm(F)的值即為蛋白質序列P和譜圖S的匹配分數.在最優路徑中，Di′j′(f′)的前驅點標記為Di j(f)，如果2個數值對(ai,bj)和(ai′,bj′)滿足條件：

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

(1)

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

通常，如果i

(2)

Di j(f)的計算為

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

(3)

Hi j(f)的狀態轉移方程為

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

(4)

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

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

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

(5)

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

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)

狀態轉移方程式(4)更新成

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

(7)

通過式(6)，可以利用動態規劃算法，最終求得全局譜圖網格中完整蛋白質序列P與譜圖S的匹配分數.

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

2 CUDA-TP算法

如第1節所述的算法雖然能在線性時間內計算出匹配分數，但當蛋白質數據庫變大或者譜圖數量增多時，其時間開銷依然不容樂觀.為了加速計算過程，本文提出了基于GPU的CUDA-TP算法，該算法首先在CPU端使用AVL樹優化前驅節點的求取，然后設計并行模型加速式(6)的計算.同時，為了加快蛋白質過濾，CUDA-TP還對蛋白質過濾算法MS-Filter[31]做了改進.算法的整體流程如圖3所示.

2.1 計算diag

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

Fig. 4 Calculate diag of all nodes in the spectra grid圖4 計算譜圖網格中所有點的diag

由式(1)易知，如果2個點是余對角的，那么它們所在坐標對應的數值對(ai,bj)的差值d在誤差范圍內是相等的.也就是說，集合{A0,A1,A2}對應相同的d.利用這條性質，可以通過AVL樹只遍歷一遍譜圖網格求取所有點的diag，如圖4所示.計算分步驟：

1) 將譜圖網格從左到右按列進行劃分，每列元素的順序從上到下；

2) 建立一棵空的AVL樹，樹的節點包含1個坐標和1個d值；

3) 從第1列的第1個元素開始，依次計算坐標元素的d值.如果AVL樹中沒有這個元素的d值，則把此值和它的坐標作為新節點插入AVL樹.如果查找到已經有節點的d值與它相等，則把該元素的diag置為此節點所存儲的坐標，同時更新節點的坐標為該元素的坐標.

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

算法1. 在CPU端計算diag.

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

輸出：所有點的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樹在查找d值時，如果AVL樹節點中的值與該值的差的絕對值小于β，就認為它們是相等的.avl_tree[d](i,j)表示AVL樹中查找到的等于d值的節點所包含的坐標.

2.2 計算diag(δ)

Fig. 5 Calculate diag(δ) of all nodes in the spectra grid圖5 計算譜圖網格中所有點的diag(δ)

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

由diag(δ)的計算過程可知，網格中的每個元素在求解時，其結果互不影響，即所有元素可以同時計算.這種性質很好地符合了CUDA并行架構要求.在CUDA架構中，CPU作為主機(host)，GPU作為主機的協處理器(co-processor)，它被看作執行高度線程化并行處理任務的計算設備(device).運行在GPU上的CUDA并行計算函數稱為內核(kernel)，以線程柵格(grid)形式組織，線程柵格由若干個線程塊(block)組成，每個線程塊又包含若干個線程(threads).實際運行時，線程塊被分割成更小的線程束(wrap)來執行運算任務，1個線程束由連續的32個線程組成.

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

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

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

輸出：存儲所有點的diag(δ)數組g.

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

②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值；行③④開始遍歷坐標左上部區域的所有點；行⑤～求當前遍歷點的delta值，如果PTM集合ΔPTM包含此值，則將其更新到數組g，最后返回的數組g即為最終結果；代碼的行⑥檢測ΔPTM中是否存在δ與delta之差的絕對值小于β，如果存在，則判斷集合ΔPTM包含delta.

2.3 數組計算并行架構

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

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

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

Fig. 7 Parallel processing architecture圖7 并行計算架構

算法3. 在GPU端計算數組E.

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

輸出：所求數組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(δ))；*獲取結果并填充進數組E*

⑧ end if

⑨ end for

⑩ end for

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

2.4 時間復雜度分析

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

(8)

因此，質譜匹配的總時間復雜度由O(N2+FN)變為

(9)

2.5 蛋白質過濾

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

1) 在低精度下計算譜圖與蛋白質的卷積數組；

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

3) 在高精度下依次對合并后的候選區域計算卷積數組，如果卷積數組中對角線分值大于設定的閾值，則將閾值更新為此分值，最后的閾值即為蛋白質的過濾分數.

蛋白質的過濾加速體現在步驟2的候選區域合并，通過合并分散的小區域來減少步驟3高精度卷積數組的個數.

CUDA-TP創建L個流(stream)對過濾出的候選蛋白質與譜圖的匹配分數并行計算，CUDA架構的流并行是粗粒度的并行，它能夠使GPU運算的同時與CPU進行數據交換.如圖8所示，為每個候選蛋白質分配1個流，流與流之間并行執行.

Fig. 8 Stream parallelism for calculating the matching score between the spectrum and candidate proteins圖8 譜圖與候選蛋白質匹配分數計算的流并行

3 實驗與結果

3.1 數據集合

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

由于TD方法是對完整蛋白質的鑒定，譜峰數量太少會導致匹配分數過低，所以只保留前體質量(precursor mass)大于2500Da且至少包含10個譜峰的譜圖[14].最終的實驗數據集及其分布如表1和圖9所示.

Table1TheNumberofSpectrafromFourDatasetsandtheNumberofProteinsintheCorrespondingDatabases

表1 數據集中的譜圖數量及數據庫中蛋白質個數

Fig. 9 Distribution of the number of peaks and the length of proteins圖9 譜峰數量和蛋白質長度分布

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

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

3.2 CUDA-TP運行時間

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

Table 2 The Running Environment表2 實驗環境

Fig. 10 Running time of CUDA-TP and MS-TopDown on four datasets with different parameters 圖10 4個數據集上CUDA-TP與MS-TopDown在不同參數下的運行時間

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

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

MS-Align+只可以調整F的大小，L默認固定為20，因此本文對比分析了3種方法在不同F參數下的運行時間，實驗結果如圖11所示:

Fig. 11 Running time of three methods on four datasets with different parameters圖11 4個數據集上3種方法在不同參數下的運行時間

3種方法的平均運行時間在表3中給出，從表3中可以看到CUDA-TP的運行速度約是MS-Align+的2倍，這證明基于譜圖比對思想的并行計算方法明顯優于通過減少譜圖網格搜索空間來加速計算過程的策略.

Table 3 Average Running Time of Three Methods表3 3種方法平均運行時間 min

3.3 CUDA-TP算法吞吐率

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

CUDA-TP的時間復雜度已經在2.4節詳細給出，3種方法在不同數據集上的吞吐率如表4所示，吞吐率指單位時間內算法執行的計算單元數量.從表4中可以看出CUDA-TP吞吐率約為MS-TopDown和MS-Align+的11倍，這在Ecoli和H4數據上表現尤為明顯，說明譜圖網格元素的增多雖然增加了運行時間，但吞吐率并沒有隨著降低.MS-Align+通過減少求取網格的數量來減少計算時間，但吞吐率并沒有實質提升.因此，CUDA-TP在運行時間和算法的吞吐率上要明顯優于目前的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在數據集上的執行時間通常要花費大量時間(本文中MS-TopDown在最大數據集上的運行時間在2 d左右)，這大大降低了其在蛋白質鑒定中的實用性，而本文提出的基于GPU的蛋白質鑒定并行算法運行速度是其10倍，有效地解決了譜圖比對思想應用于大規模數據的弊端，與現有的MS-Align+算法相比具有明顯優勢，通過上述多種不同的實驗測試，驗證了CUDA-TP的優秀性能.CUDA-TP源代碼托管在github公共服務網站①.

4 總結及展望

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

在TD策略下對完整蛋白質進行鑒定，除了計算蛋白質與譜圖的匹配分數，還需要進一步評估匹配結果的統計顯著性.因此如何計算匹配分數的同時得到蛋白質與譜圖匹配的統計顯著性評估結果，并且不降低時間運行效率，這將是我們下一步的研究工作.