鐘誠(chéng)，孫輝

（1.廣西大學(xué)計(jì)算機(jī)與電子信息學(xué)院，廣西 南寧 530004；2.廣西高校并行分布與智能計(jì)算重點(diǎn)實(shí)驗(yàn)室，廣西 南寧 530004）

0 引言

基因組數(shù)據(jù)是一類(lèi)重要的醫(yī)療大數(shù)據(jù)，被廣泛應(yīng)用于生物醫(yī)學(xué)、定制化醫(yī)療服務(wù)、法醫(yī)學(xué)鑒定和全基因組關(guān)聯(lián)研究等領(lǐng)域[1-2]。基因組數(shù)據(jù)敏感序列是指可用于發(fā)起推斷攻擊、成員身份攻擊、身份追蹤攻擊和完全攻擊等基因組數(shù)據(jù)隱私攻擊的核苷酸序列，包括短串聯(lián)重復(fù)（STR,short tandem repeat）序列和疾病相關(guān)序列（DRS,disease-related sequence）[3]。不同于網(wǎng)絡(luò)隱私數(shù)據(jù)的可重塑性，基因組數(shù)據(jù)敏感序列信息的泄露會(huì)對(duì)隱私泄露者的醫(yī)療保險(xiǎn)、社會(huì)就業(yè)和公共福利等產(chǎn)生不利影響[4]。通過(guò)識(shí)別過(guò)濾基因組數(shù)據(jù)中的敏感信息，可以實(shí)現(xiàn)基因組數(shù)據(jù)敏感信息脫敏，有助于降低基因組數(shù)據(jù)隱私敏感信息受攻擊的風(fēng)險(xiǎn)。

基因組數(shù)據(jù)中的敏感序列包括STR 和DRS 這2 類(lèi)。敏感序列識(shí)別串行算法的設(shè)計(jì)主要基于Hash映射思想，即通過(guò)產(chǎn)生序列Hash 指印來(lái)實(shí)現(xiàn)敏感序列的精確識(shí)別。此類(lèi)代表算法有SRF（short read filter）[3]、LRF（long read filter）[4]和PriLive[5]等。這些基于序列Hash 指印思想的串行算法對(duì)高錯(cuò)誤率的長(zhǎng)序列基因組數(shù)據(jù)敏感序列識(shí)別準(zhǔn)確率較低，且難以滿(mǎn)足對(duì)大規(guī)模長(zhǎng)序列實(shí)時(shí)識(shí)別過(guò)濾的要求。為此，一些學(xué)者基于多核CPU 集群系統(tǒng)，通過(guò)設(shè)計(jì)多核CPU 并行算法，在準(zhǔn)確識(shí)別含有錯(cuò)誤信息的敏感序列的同時(shí)，加速求解大規(guī)?；蚪M數(shù)據(jù)敏感序列識(shí)別問(wèn)題。

關(guān)于基因組數(shù)據(jù)敏感序列識(shí)別并行算法的研究，文獻(xiàn)[6]基于STR 挖掘工具FLASH（fast length adjustment of short read）和PERF（Python exhaustive repeat finder），在Hadoop 平臺(tái)上開(kāi)發(fā)了CPU 并行加速的基因組數(shù)據(jù)短串聯(lián)重復(fù)序列識(shí)別算法BigFiRSt。文獻(xiàn)[7]采用劃分后綴分組方法，提出基于CPU 集群系統(tǒng)的最大串聯(lián)重復(fù)序列識(shí)別并行算法。在串行算法FindTRs 的基礎(chǔ)上，文獻(xiàn)[8]設(shè)計(jì)了在多核CPU 集群系統(tǒng)上運(yùn)行的大規(guī)模短串聯(lián)重復(fù)序列識(shí)別并行算法MPI-dot2dot。根據(jù)STR 的周期重復(fù)性質(zhì)，文獻(xiàn)[9]在多核CPU 集群系統(tǒng)上設(shè)計(jì)實(shí)現(xiàn)了生物序列模式精確匹配的并行算法。文獻(xiàn)[10]通過(guò)枚舉所有具有統(tǒng)計(jì)顯著性的局部最優(yōu)出現(xiàn)頻率的子串，提出了一種面向多核CPU 系統(tǒng)的頻繁子串模式挖掘并行算法。上述多核CPU 并行算法僅針對(duì)基因組數(shù)據(jù)中的STR 進(jìn)行識(shí)別，并沒(méi)有處理基因組數(shù)據(jù)中DRS 的識(shí)別，因此這些并行算法僅能識(shí)別出基因組數(shù)據(jù)中的部分敏感序列。對(duì)于低錯(cuò)誤率的基因組數(shù)據(jù)，文獻(xiàn)[11]提出一種多核CPU 并行算法，以識(shí)別基因組數(shù)據(jù)中的短串聯(lián)重復(fù)序列和疾病相關(guān)序列。

已有的基因組數(shù)據(jù)敏感序列識(shí)別并行算法的設(shè)計(jì)主要是針對(duì)敏感序列2 種類(lèi)型（STR 和DRS）中的一種進(jìn)行處理的；有的并行算法雖然同時(shí)處理STR 和DRS 的識(shí)別，但它面向的是低錯(cuò)誤率序列構(gòu)成的基因組數(shù)據(jù)，而不是高錯(cuò)誤率序列構(gòu)成的基因組數(shù)據(jù)。此外，已有的基因組數(shù)據(jù)敏感序列識(shí)別并行算法是多核CPU 并行算法，隨著單分子實(shí)時(shí)測(cè)序和牛津納米孔測(cè)序等第三代測(cè)序技術(shù)的發(fā)展，測(cè)序平臺(tái)產(chǎn)生的序列具有長(zhǎng)度長(zhǎng)和錯(cuò)誤率高的特點(diǎn)[12]。已有的并行算法難以高效、準(zhǔn)確地識(shí)別出高錯(cuò)誤率長(zhǎng)序列構(gòu)成的基因組數(shù)據(jù)中的2 類(lèi)敏感序列STR 和DRS。為解決此問(wèn)題，本文研究改進(jìn)序列相似度計(jì)算模型，進(jìn)而提出一種CPU 和GPU 協(xié)同計(jì)算加速識(shí)別高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)中STR 和DRS 的并行算法CGPU-F3SR。該并行算法的特色和創(chuàng)新介紹如下。

1) 采用滑動(dòng)窗口策略和k-mer 編碼方法，分割并提取高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)中的核苷酸序列堿基錯(cuò)誤信息。引入布隆過(guò)濾器，對(duì)高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)進(jìn)行高效過(guò)濾，避免算法對(duì)重復(fù)數(shù)據(jù)的多次計(jì)算。

2) 根據(jù)短串聯(lián)重復(fù)序列和疾病相關(guān)序列的結(jié)構(gòu)特點(diǎn)，提出一種改進(jìn)的序列相似度計(jì)算模型，實(shí)現(xiàn)準(zhǔn)確識(shí)別出高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)中的2 類(lèi)敏感序列（STR 和DRS）。

3) 對(duì)基因組數(shù)據(jù)測(cè)序序列進(jìn)行周期劃分，設(shè)計(jì)實(shí)現(xiàn)CPU 和GPU 協(xié)同計(jì)算的敏感序列識(shí)別并行算法，大大提升算法的運(yùn)行效率。

1 方法

設(shè)待識(shí)別基因組數(shù)據(jù)共有l(wèi)條長(zhǎng)序列，Ns[0:ns-1]表示長(zhǎng)度為ns的第s條長(zhǎng)序列，s=0,1,2,…,l-1。本文提出的CPU 和GPU 協(xié)同計(jì)算的高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)敏感序列識(shí)別并行算法包括序列分割過(guò)濾、短串聯(lián)重復(fù)序列并行識(shí)別、疾病相關(guān)序列并行識(shí)別和生成掩碼核苷酸序列4 個(gè)階段。

1.1 并行算法4 個(gè)階段

1.1.1 序列分割過(guò)濾

為并行識(shí)別高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)中的敏感序列，首先將長(zhǎng)序列Ns[0:ns-1]分割為ns-r+1條、每條長(zhǎng)度為r的短序列Rs,j=Ns[j:j+r-1]，其中，短序列Rs,j和Rs,j+1共享r-1 個(gè)重疊堿基，j=0,1,2,…,ns-r，s=0,1,2,…,l-1，ns＞r。因?yàn)槿祟?lèi)DNA序列中約99.5%的序列具有高度一致性，所以將長(zhǎng)序列分割為短序列會(huì)產(chǎn)生許多相同的短序列[3]。

布隆過(guò)濾器由一個(gè)二進(jìn)制位數(shù)組（位向量）和若干Hash 函數(shù)組成，多用于數(shù)據(jù)去重以及數(shù)據(jù)檢索[13]。為避免對(duì)重復(fù)短序列的多次計(jì)算，本文引入雙布隆過(guò)濾機(jī)制對(duì)短序列Rs,j并行過(guò)濾去重，j=0,1,2,…,ns-r，s=0,1,2,…,l-1。

設(shè)布隆過(guò)濾器BF1用于敏感序列去重，其位數(shù)組和使用的Hash 函數(shù)數(shù)目分別為B1[0:b1-1]和h1。布隆過(guò)濾器BF2用于非敏感序列去重，其位數(shù)組和使用的Hash 函數(shù)數(shù)目分別為B2[0:b2-1]和h2。為了初始化雙布隆過(guò)濾器，需根據(jù)式(1)計(jì)算B1和B2的位數(shù)組規(guī)模b1和b2[13]。

其中，pi表示布隆過(guò)濾器誤報(bào)率，ti表示BFi可容納序列數(shù)目的上限。通過(guò)預(yù)設(shè)參數(shù)pi和ti，可以根據(jù)式(1)計(jì)算得到位數(shù)組規(guī)模bi，i=1,2。此時(shí)，可以計(jì)算BFi使用的Hash 函數(shù)數(shù)目hi[13]為

依據(jù)式(1)和式(2)計(jì)算出雙布隆過(guò)濾器的參數(shù)b1和b2以及h1和h2之后，初始化布隆過(guò)濾器位數(shù)組B1,i=0 和B2,k=0，i=0,1,2,…,b1-1，k=0,1,2,…,b2-1。

當(dāng)初始化雙布隆過(guò)濾器完成后，將短序列Rs,j進(jìn)行周期劃分，然后使用Pr 個(gè)CPU 核心通過(guò)雙布隆過(guò)濾器BF1和BF2并行地對(duì)短序列Rs,j進(jìn)行過(guò)濾去重，s=0,1,2,…,l-1，j=0,1,2,…,ns-r。

為方便理解短序列周期劃分思想，圖1 給出了塊劃分和周期劃分短序列的示例。

人類(lèi)基因組數(shù)據(jù)測(cè)序序列中存在連續(xù)的串聯(lián)重復(fù)區(qū)域。若將分割得到的短序列采用圖1(a)所示塊劃分方式并行處理，則會(huì)造成部分線程處于“計(jì)算等待”狀態(tài)。為均衡數(shù)據(jù)計(jì)算，使每個(gè)線程達(dá)到負(fù)載均衡，本文采用圖1(b)所示的周期劃分對(duì)短序列進(jìn)行并行處理，并采用雙布隆過(guò)濾器過(guò)濾去重。

圖1 基于塊劃分和基于周期劃分短序列的示例

本文并行短序列過(guò)濾去重方法介紹如下。

Pr 個(gè)CPU 核心并行對(duì)Rs,v+dPr執(zhí)行h1次Hash計(jì)算，以獲得h1個(gè)Hash 值H1,i(Rs,v+dPr)，且對(duì)BF1中的位數(shù)組B1進(jìn)行查詢(xún)。

若B1[H1,i(Rs,v+dPr)]的值均為1，i=0,1,2,…,h1-1，則表明Rs,v+dPr是重復(fù)短序列且為敏感短序列；若B1[H1,i(Rs,v+dPr)]的值不全為1，則需對(duì)Rs,v+dPr執(zhí)行h2次Hash 計(jì)算，以獲得h2個(gè)Hash 值H2,i(Rs,v+dPr)，并對(duì)BF2中的位數(shù)組B2進(jìn)行查詢(xún)。

若B2[H2,i(Rs,v+dPr)]的值均為1，i=0,1,2,…,h2-1，則表明Rs,v+dPr為重復(fù)短序列且為非敏感短序列；若B2[H2,i(Rs,v+dPr)]的值不全為1，則表明Rs,v+dPr不存在于BF1或BF2中，Rs,v+dPr不是一條重復(fù)序列，此時(shí)需對(duì)Rs,v+dPr進(jìn)行短串聯(lián)重復(fù)和疾病相關(guān)序列識(shí)別處理，d=0,1,2,…,，v+dPr≤ns-r，v=0,1,2,…,Pr-1，s=0,1,2,…,l-1。

通過(guò)并行地對(duì)短序列過(guò)濾去重，可避免對(duì)短序列的重復(fù)計(jì)算。若短序列不存在于布隆過(guò)濾器BF1或BF2，則需對(duì)短序列進(jìn)行短串聯(lián)重復(fù)序列識(shí)別的處理。

1.1.2 短串聯(lián)重復(fù)序列并行識(shí)別

STR 也被稱(chēng)作微衛(wèi)星DNA，由1～6 個(gè)堿基的重復(fù)基序組成，約占人類(lèi)基因組數(shù)據(jù)的3%[3]。短串聯(lián)重復(fù)序列的重復(fù)結(jié)構(gòu)性質(zhì)表明，通過(guò)合理分割短串聯(lián)重復(fù)序列，可以得到記錄短串聯(lián)重復(fù)序列基本重復(fù)單元的相同子段[14]。在借鑒文獻(xiàn)[14]思想的基礎(chǔ)上，本文使用Pr 個(gè)CPU 核心并行地識(shí)別高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)中的短串聯(lián)重復(fù)序列。

并行識(shí)別短串聯(lián)重復(fù)序列的方法介紹如下。

首先，Pr 個(gè)CPU 核心并行分割短序列Rs,v+dPr為4 個(gè)長(zhǎng)度均為w的短片段集合i=0,1,2,3，v=0,1,2,…,Pr-1，v+dPr≤ns-r，s=0,1,…,l-1。

為降低將非敏感序列誤識(shí)別為短串聯(lián)重復(fù)序列的概率，根據(jù)短串聯(lián)重復(fù)序列的結(jié)構(gòu)特點(diǎn)，文獻(xiàn)[15]引入控制參數(shù)c_str 和皮爾遜相關(guān)系數(shù)（PCC,Pearson correlation coefficient）改進(jìn)序列局部片段相似度的計(jì)算為

1.1.3 疾病相關(guān)序列并行識(shí)別

DRS 是指攜帶疾病易感基因的核苷酸序列[11]。不同于短串聯(lián)重復(fù)序列，疾病相關(guān)序列更多體現(xiàn)在位點(diǎn)上的變異，而非結(jié)構(gòu)上的重復(fù)。為此，本文通過(guò)改進(jìn)序列相似度計(jì)算模型，將待識(shí)別序列和第三方數(shù)據(jù)庫(kù)中疾病相關(guān)序列進(jìn)行計(jì)算比對(duì)，以并行識(shí)別出高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)中的疾病相關(guān)序列。

1) 并行構(gòu)建敏感序列詞典

構(gòu)建敏感序列詞典是為了降低大量短序列與第三方數(shù)據(jù)庫(kù)中疾病相關(guān)序列相似度計(jì)算的代價(jià)。本文使用二維數(shù)組 hostDict[0:m-1][0:4k+2]表示CPU 端構(gòu)建包含m條序列的敏感序列詞典。

并行構(gòu)建敏感序列詞典的方法介紹如下。

首先，Pr 個(gè)CPU 核心并行地從第三方數(shù)據(jù)庫(kù)提取每個(gè)疾病易感基因的染色體編號(hào)CHR_ID 和染色體位置CHR_POS，并使用pysam 工具[16]從人類(lèi)參考基因組HG38 中提取m條長(zhǎng)度為r的疾病相關(guān)序列Dv,i，i=0,1,2,…,m-1，v=0,1,2,…,Pr-1。

2) 疾病相關(guān)序列相似度并行計(jì)算

在判斷短序列Rs,j是否為疾病相關(guān)序列之前，本文已經(jīng)對(duì)Rs,j進(jìn)行過(guò)濾去重和短串聯(lián)重復(fù)序列識(shí)別。因此，僅需要對(duì)經(jīng)上述步驟處理后剩余的g條序列進(jìn)行疾病相關(guān)序列識(shí)別。將剩余的g條序列和敏感序列詞典hostDict[0:m-1][0:4k+2]中的m條序列進(jìn)行相似度計(jì)算比對(duì)十分耗時(shí)。為此，本文采用多核CPU 和GPU 協(xié)同并行計(jì)算以加速疾病相關(guān)序列的識(shí)別，j=0,1,2,…,ns-r，s=0,1,2,…,l-1。

CUDA（computer unified device architecture）并行編程模型將GPU 線程層次結(jié)構(gòu)抽象成線程塊網(wǎng)格Grid 和線程塊Block 兩層[17]。在CUDA 線程層次結(jié)構(gòu)中，Block 和Grid 由uint3 類(lèi)型定義的三維向量 blockIdx.x、blockIdx.y、blockIdx.z 和threadIdx.x、threadIdx.y、threadIdx.z 表示[17]。為準(zhǔn)確定位識(shí)別疾病相關(guān)序列，本文采用一維線程塊網(wǎng)格和線程塊進(jìn)行CUDA 線程索引。

其中，e=0,1,2，i=0,1,…,m-1，tID=0,1,…,g-1，s=0,1,…,l-1，dis_sim 為疾病相關(guān)序列相似度閾值。若sim(Di,Rs,tID)≥dis_sim，則序列Rs,tID是一條疾病相關(guān)序列，且它和敏感序列詞典中第i條序列Di具有很高的序列相似度。

1.1.4 生成掩碼核苷酸序列

1.2 并行算法描述與分析

設(shè)GWAS Catlog 數(shù)據(jù)庫(kù)中疾病相關(guān)序列數(shù)據(jù)集為disData（包含m組數(shù)據(jù)），HG 表示疾病相關(guān)序列并行識(shí)別步驟構(gòu)建敏感序列詞典時(shí)所使用的參考基因組，并行算法使用的CPU 核心數(shù)為Pr，使用的GPU 線程數(shù)為Cu。

圖2 給出了CPU/GPU 并行算法CGPU-F3SR的處理流程。

圖2 中，res[0:ns-1]表示記錄每條短序列的敏感序列識(shí)別結(jié)果的輔助數(shù)組，且res[j]∈{0,1,2,3,4}分別表示原始短序列結(jié)果標(biāo)記、BF1過(guò)濾標(biāo)記、BF2過(guò)濾標(biāo)記、短串聯(lián)重復(fù)序列標(biāo)記和疾病相關(guān)序列標(biāo)記，j=0,1,2,…,ns-1，s=0,1,2,…,l-1。hostRes[0:g-1]表示保存疾病相關(guān)序列的識(shí)別結(jié)果的輔助數(shù)組，其中，，I(res[j])=0 表示指示函數(shù)（條件成立則取值為1，否則取值為0）。strIndex[0:g-1]表示保存疾病相關(guān)序列識(shí)別的短序列索引的輔助數(shù)組。式(4)中，的求解需要做許多分支判斷。若使用GPU 計(jì)算，則將使部分GPU 線程處于長(zhǎng)時(shí)間“計(jì)算等待”狀態(tài)，為均衡 GPU 端數(shù)據(jù)運(yùn)算，本文引入輔助數(shù)組hostData[0:g-1][0:4k+2]，并使用多核CPU 預(yù)編碼核苷酸序列和堿基含量統(tǒng)計(jì)，以簡(jiǎn)化GPU 端數(shù)據(jù)運(yùn)算，提升GPU 端的運(yùn)行效率。

圖2 CPU/GPU 并行算法CGPU-F3SR 的處理流程

算法1 給出了CGPU-F3SR 的流程。

算法1 CGPU-F3SR

當(dāng)CPU/GPU 并行算法CGPU-F3SR 中參數(shù)的gpuBlock 取值為0 時(shí)，它就變?yōu)閮H采用多核CPU并行計(jì)算的版本（簡(jiǎn)記為算法CPU-F3SR）。

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

實(shí)驗(yàn)在廣西大學(xué)的Sugon 7000A 超級(jí)并行計(jì)算機(jī)上進(jìn)行，使用CPU/GPU 異構(gòu)計(jì)算節(jié)點(diǎn)。每個(gè)計(jì)算節(jié)點(diǎn)的內(nèi)存容量為512 GB、外部存儲(chǔ)容量為8×900 GB、CPU 為2×Intel Xeon Gold 6230（總內(nèi)核數(shù)為40）、GPU 為NVIDIA Tesla-T4（顯存容量為16 GB，含有2 560 個(gè)CUDA 核心）。運(yùn)行操作系統(tǒng)為64 位版本CentOS7.4。采用C++、CUDA 和OpenMP 混合編程實(shí)現(xiàn)并行算法CGPU-F3SR。

實(shí)驗(yàn)采用GWAS Catalog 數(shù)據(jù)庫(kù)[18]中全基因組關(guān)聯(lián)研究數(shù)據(jù)集gwas_catalog_v1.0.tsv構(gòu)建敏感序列詞典，該數(shù)據(jù)集包含Scoliosis、HIV-1 replication 和Opioid sensitivity 等4 680 種疾病，共計(jì)216 250 組數(shù)據(jù)。采用的參考基因組為NCBI 數(shù)據(jù)庫(kù)[19]開(kāi)放的人類(lèi)基因組HG38。實(shí)驗(yàn)所用的疾病易感基因來(lái)自GWAS Catalog 數(shù)據(jù)庫(kù)[18]，短串聯(lián)重復(fù)序列來(lái)自TRDB 數(shù)據(jù)庫(kù)[20]。參考文獻(xiàn)[3]的數(shù)據(jù)預(yù)處理規(guī)則，預(yù)處理每條序列長(zhǎng)度為50 bp 的基準(zhǔn)數(shù)據(jù)序列集All-Together（2 239 340 條序列），每個(gè)基準(zhǔn)數(shù)據(jù)集的非敏感序列和敏感序列數(shù)量之比為7:1。

為評(píng)估算法識(shí)別效果，需獲取不同長(zhǎng)度且?guī)в忻舾行蛄袠?biāo)簽的測(cè)序長(zhǎng)序列。為此，實(shí)驗(yàn)對(duì)錯(cuò)誤率為2%～20%的All-Together 基準(zhǔn)數(shù)序列集進(jìn)行隨機(jī)抽樣，生成序列長(zhǎng)度為100～400 kbp 的總規(guī)模約24 GB 的4 組長(zhǎng)序列數(shù)據(jù)集DataSet-100 kbp～Data Set-400 kbp 用于并行算法參數(shù)調(diào)優(yōu)測(cè)試。

實(shí)驗(yàn)將本文CPU 和GPU 協(xié)同計(jì)算的并行算法CGPU-F3SR與并行算法PARA-LRF[11]、僅采用CPU并行計(jì)算的版本CPU-F3SR（gpuBlock=0）、串行算法F3SR[15]進(jìn)行測(cè)試比較識(shí)別效果。實(shí)驗(yàn)評(píng)估算法的識(shí)別準(zhǔn)確率（acc,accuracy）[3]、查準(zhǔn)率（pre,precision）[4]和假陽(yáng)性率（fpr,false positive rate）[4]。準(zhǔn)確率和查準(zhǔn)率越高，算法識(shí)別效果越好；假陽(yáng)性率越低，算法識(shí)別效果越好。在運(yùn)行效率方面，實(shí)驗(yàn)評(píng)估了算法的識(shí)別吞吐量（thp,throughput）[4]。吞吐量定義為算法單位時(shí)間識(shí)別脫氧核糖核苷酸的數(shù)量，單位為nt/s[3-4]，吞吐量越高，算法運(yùn)行效率越高。

下面，首先給出本文所提出的CPU/GPU 并行算法CGPU-F3SR 的參數(shù)調(diào)優(yōu)實(shí)驗(yàn)結(jié)果，然后給出消融實(shí)驗(yàn)結(jié)果，最后給出算法CGPU-F3SR 和同類(lèi)并行算法的實(shí)驗(yàn)結(jié)果及分析。

2.1 并行算法參數(shù)調(diào)優(yōu)實(shí)驗(yàn)

參考文獻(xiàn)[3,4,15]給出的基因組數(shù)據(jù)敏感序列識(shí)別算法研究的相關(guān)實(shí)驗(yàn)，算法CGPU-F3SR 選取pi=0.000 1、ti=108、r=48 bp、k=5、str_sim=0.64、dis_sim=0.80 進(jìn)行算法部分參數(shù)的初始化，i=1,2。

為確定GPU 線程塊大小對(duì)并行算法吞吐量的影響，實(shí)驗(yàn)測(cè)試了并行算法CGPU-F3SR 采用不同的GPU 線程塊大小在4 組長(zhǎng)序列數(shù)據(jù)集上運(yùn)行時(shí)的吞吐量，結(jié)果如圖3 所示。

圖3 并行算法CGPU-F3SR 采用不同的GPU 線程塊大小在4 組長(zhǎng)序列數(shù)據(jù)集上運(yùn)行時(shí)的吞吐量

圖3 實(shí)驗(yàn)結(jié)果表明，當(dāng)GPU 線程塊大小固定時(shí)，序列數(shù)據(jù)集規(guī)模越大，并行算法CGPU-F3SR 的吞吐量越高。原因如下：1) 測(cè)序長(zhǎng)序列數(shù)據(jù)集規(guī)模越大，CGPU-F3SR 分割長(zhǎng)序列得到的短序列中的重復(fù)短序列越多，由于雙布隆過(guò)濾避免了重復(fù)短序列的多次運(yùn)算，因此過(guò)濾效果明顯；2) CGPU-F3SR 將編碼后的短序列提交給GPU 進(jìn)行疾病相關(guān)序列并行識(shí)別，測(cè)序數(shù)據(jù)集越大，GPU 線程平均計(jì)算處理的數(shù)據(jù)越多，算法的并行度越高。按照CUDA 并行模型，網(wǎng)格Grid中的各線程塊Block 會(huì)被分配到GPU 的各流式處理器中執(zhí)行。圖3 的實(shí)驗(yàn)結(jié)果還表明，當(dāng)參數(shù)gpuBlock取值為128～384 時(shí)，算法整體上具有較高的吞吐量。

圖4 進(jìn)一步給出了GPU 線程塊大小gpuBlock=256、使用不同 CPU 核心數(shù)運(yùn)行并行算法CGPU-F3SR 時(shí)得到的吞吐量結(jié)果。

圖4 實(shí)驗(yàn)結(jié)果表明，并行算法CGPU-F3SR 的識(shí)別吞吐量與使用的CPU 核心數(shù)呈正相關(guān)。數(shù)據(jù)集規(guī)模越大，序列越長(zhǎng)，算法吞吐量受CPU 核心數(shù)增益效果越明顯。

圖4 并行算法CGPU-F3SR 使用不同CPU 核心數(shù)在4 組長(zhǎng)序列數(shù)據(jù)集上運(yùn)行時(shí)的吞吐量

為使本文并行算法CGPU-F3SR 充分發(fā)揮CPU和GPU 協(xié)同并行計(jì)算的能力，根據(jù)圖3 和圖4 的實(shí)驗(yàn)結(jié)果，本文選取gpuBlock=256 和Pr=4 的組合進(jìn)行后續(xù)的實(shí)驗(yàn)。

2.2 并行算法消融實(shí)驗(yàn)

并行算法CGPU-F3SR 引入布隆過(guò)濾器對(duì)序列數(shù)據(jù)進(jìn)行過(guò)濾，以免對(duì)分割短序列的重復(fù)計(jì)算。為了確定過(guò)濾機(jī)制對(duì)并行算法CGPU-F3SR 的性能影響，本文實(shí)驗(yàn)評(píng)估了并行算法CGPU-F3SR 進(jìn)行序列過(guò)濾With Filter（參數(shù)Pr=4，gpuBlock=256）和不進(jìn)行序列過(guò)濾 Without Filter（參數(shù) Pr=30，gpuBlock=256）2 種情形在4 組長(zhǎng)序列數(shù)據(jù)集運(yùn)行得到的識(shí)別準(zhǔn)確率、查準(zhǔn)率、假陽(yáng)性率和吞吐量，實(shí)驗(yàn)結(jié)果如表1 所示。

表1 的消融實(shí)驗(yàn)結(jié)果表明，在識(shí)別性能方面，并行算法CGPU-F3SR 在With Filter 和Without Filter 這2 種情形下均獲得相同的識(shí)別準(zhǔn)確率、查準(zhǔn)率和假陽(yáng)性率。這是因?yàn)椋?) 布隆過(guò)濾器的多哈希映射機(jī)制可以對(duì)相同序列產(chǎn)生多個(gè)相同的哈希映射，這可以在避免對(duì)重復(fù)序列的多次計(jì)算的同時(shí)不影響其識(shí)別性能；2) 改進(jìn)序列相似度計(jì)算模型能對(duì)含有錯(cuò)誤信息的相似序列產(chǎn)生近似度量值，并行算法CGPU-F3SR 的識(shí)別準(zhǔn)確率、查準(zhǔn)率和假陽(yáng)性率均取決于改進(jìn)的序列相似度計(jì)算模型，而不是取決于布隆過(guò)濾器。在運(yùn)行效率方面，算法使用4 個(gè)CPU核心的With Filter 模式的識(shí)別吞吐量比使用30 個(gè)CPU 核心的Without Filter 模式高了2 個(gè)數(shù)量級(jí)。這是因?yàn)殚L(zhǎng)序列基因組數(shù)據(jù)中分割短序列具有大量相同的重復(fù)短序列，布隆過(guò)濾效果明顯。

表1 并行算法CGPU-F3SR 的識(shí)別準(zhǔn)確率、查準(zhǔn)率、假陽(yáng)性率和吞吐量

2.3 并行算法識(shí)別性能實(shí)驗(yàn)

為評(píng)估算法的并行化對(duì)高錯(cuò)誤率長(zhǎng)序列基因數(shù)據(jù)中敏感序列（包括STR 和DRS）識(shí)別效果的影響，首先，實(shí)驗(yàn)測(cè)試CPU 和GPU 協(xié)同并行算法CGPU-F3SR、僅CPU 并行的算法CPU-F3SR 和串行算法F3SR（參數(shù)取值pi=0.000 1，ti=104，r=48，k=5，dis_sim=0.80，str_sim=0.64，i=1,2）[15]在4 組長(zhǎng)序列數(shù)據(jù)集上運(yùn)行得到的識(shí)別準(zhǔn)確率、查準(zhǔn)率和假陽(yáng)性率，實(shí)驗(yàn)結(jié)果如表2 所示。

表2 的實(shí)驗(yàn)結(jié)果表明，CPU 與GPU 協(xié)同并行算法CGPU-F3SR 和CPU 并行化算法CPU-F3SR 均獲得了與串行算法F3SR 相同的識(shí)別結(jié)果。這是因?yàn)殚L(zhǎng)序列分割得到的相鄰短序列之間共享r-1 個(gè)高度重疊的堿基，通過(guò)CPU/GPU 并行計(jì)算不會(huì)影響算法的最終識(shí)別質(zhì)量。在上述實(shí)驗(yàn)的基礎(chǔ)上，進(jìn)一步測(cè)試比較本文并行算法CGPU-F3SR 和同類(lèi)并行算法PARA-LRF（參數(shù)p=0.02，m=68 915 861）[11]在DataSet-100 kbp～DataSet-400 kbp 數(shù)據(jù)集上運(yùn)行時(shí)，對(duì)高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)中敏感序列（包括STR 和DRS）進(jìn)行識(shí)別獲得的識(shí)別準(zhǔn)確率、查準(zhǔn)率和假陽(yáng)性率，實(shí)驗(yàn)結(jié)果如表3 所示。

表2 并行算法CGPU-F3SR、CPU-F3SR 和F3SR 的識(shí)別準(zhǔn)確率、查準(zhǔn)率和假陽(yáng)性率

表3 的實(shí)驗(yàn)結(jié)果表明，相較于并行算法PARA-LRF，本文CPU/GPU 并行算法CGPU-F3SR整體上具有更高的識(shí)別準(zhǔn)確率和查準(zhǔn)率以及更低的識(shí)別假陽(yáng)性率。這是因?yàn)椴⑿兴惴–GPU-F3SR采用k-mer 編碼策略有效地提取了高錯(cuò)誤率長(zhǎng)序列中的堿基錯(cuò)誤信息，且根據(jù)STR 和DRS 的結(jié)構(gòu)特點(diǎn)，通過(guò)引入控制參數(shù)改進(jìn)了STR 和DRS 相似度的計(jì)算方法，使本文并行算法CGPU-F3SR 在準(zhǔn)確識(shí)別出高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)中敏感序列的同時(shí)，降低了將非敏感序列誤判為敏感序列的情形的發(fā)生，從而整體上獲得更好的識(shí)別效果。

表3 并行算法CGPU-F3SR 和PARA-LRF 的識(shí)別準(zhǔn)確率、查準(zhǔn)率和假陽(yáng)性率

2.4 并行算法運(yùn)行效率實(shí)驗(yàn)

在基因組數(shù)據(jù)敏感序列識(shí)別研究中，算法吞吐量和運(yùn)行時(shí)間相關(guān)，多用于運(yùn)行效率評(píng)估?；蚪M數(shù)據(jù)敏感序列識(shí)別并行算法相關(guān)研究較少，目前缺乏CPU 和 GPU 協(xié)同計(jì)算的并行算法。并行算法PARA-LRF[11]是基于多核CPU 設(shè)計(jì)的并行算法。為公平起見(jiàn)，實(shí)驗(yàn)首先測(cè)試了本文僅多核CPU 并行算法CPU-F3SR和多核CPU并行算法PARA-LRF[11]使用不同CPU核心數(shù)在4組長(zhǎng)序列數(shù)據(jù)集上運(yùn)行得到的識(shí)別吞吐量，結(jié)果如圖5 所示。

圖5 多核CPU 并行算法CPU-F3SR 和PARA-LRF 使用的CPU 核心數(shù)對(duì)吞吐量的影響

從圖5(a)可以看出，在DataSet-100 kbp 數(shù)據(jù)集上，當(dāng)使用的CPU 核心數(shù)為1～28 時(shí)，CPU 并行算法PARA-LRF 的識(shí)別吞吐量高于本文CPU 并行算法CPU-F3SR；當(dāng)使用的CPU 核心數(shù)為28～40時(shí)，并行算法CPU-F3SR 比PARA-LRF 獲得了整體上更高的識(shí)別吞吐量。

從圖5(b)～圖5(d)可以看出，對(duì)序列長(zhǎng)度為200～400 kbp 的3 個(gè)長(zhǎng)序列數(shù)據(jù)集，當(dāng)CPU 核心數(shù)分別為22、24 和26 時(shí)，并行算法CPU-F3SR 達(dá)到吞吐量拐點(diǎn)，此時(shí)本文多核 CPU 并行算法CPU-F3SR 比多核CPU 并行算法PARA-LRF 整體上具有更高的識(shí)別吞吐量。這是因?yàn)椴⑿兴惴–PU-F3SR 并行過(guò)濾去重分割得到的短序列、并行識(shí)別短串聯(lián)重復(fù)序列、并行識(shí)別疾病相關(guān)序列的部分占比高、算法并行度高，所以在4 組長(zhǎng)序列數(shù)據(jù)集上運(yùn)行時(shí)獲得了更高的識(shí)別吞吐量。

從長(zhǎng)序列基因組數(shù)據(jù)中識(shí)別疾病相關(guān)序列的過(guò)程是一個(gè)計(jì)算密集型工作，為進(jìn)一步加速求解敏感序列識(shí)別問(wèn)題，本文設(shè)計(jì)實(shí)現(xiàn)了CPU 和GPU 協(xié)同計(jì)算的敏感序列識(shí)別并行算法CGPU-F3SR。為此，實(shí)驗(yàn)也測(cè)試了 CPU/GPU 并行算法CGPU-F3SR（參數(shù)Pr=4、gpuBlock=256）、多核CPU并行算法PARA-LRF（參數(shù)p=0.02、m=68 915 861、Pr=4）[11]、多核CPU 并行算法CPU-F3SR（參數(shù)Pr=4）在4 組長(zhǎng)序列數(shù)據(jù)集上運(yùn)行得到的識(shí)別吞吐量，結(jié)果如表4 所示。

從表4 可以看出，相較于本文僅采用多核CPU并行加速的CPU-F3SR 算法，本文的CPU 和GPU協(xié)同并行求解敏感序列識(shí)別問(wèn)題的CGPU-F3SR 算法平均吞吐量高達(dá)11.033×106nt/s，極大提升了識(shí)別吞吐量；相較于已有的多核 CPU 并行算法PARA-LRF，本文算法CGPU-F3SR 在吞吐量方面也具有較大幅度的提升。CGPU-F3SR 算法的平均吞吐量分別是CPU-F3SR 和PARA-LRF 算法的11.45 倍和2.44 倍。CGPU-F3SR 算法之所以獲得更高的吞吐量，是因?yàn)樗捎肎PU 并行加速處理屬于計(jì)算密集型的疾病相關(guān)序列識(shí)別過(guò)程，加速效果顯著。

表4 并行算法CGPU-F3SR、CPU-F3SR 和PARA-LRF 的識(shí)別吞吐量

綜上所述，CPU 和GPU 協(xié)同計(jì)算的并行算法CGPU-F3SR 在準(zhǔn)確識(shí)別出大規(guī)模高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)中的2 類(lèi)敏感序列（STR 和DRS）的同時(shí)，顯著提升了吞吐量。

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

第三代測(cè)序技術(shù)產(chǎn)生的基因組測(cè)序數(shù)據(jù)正以超摩爾定律的速度顯著增長(zhǎng)。隨著人們生物信息安全意識(shí)的逐漸增強(qiáng)，保護(hù)個(gè)體基因組數(shù)據(jù)敏感序列變得越來(lái)越重要。本文提出的CPU 和GPU 協(xié)同計(jì)算識(shí)別的并行算法CGPU-F3SR 的特色和優(yōu)勢(shì)是高效準(zhǔn)確地識(shí)別出第三代測(cè)序技術(shù)產(chǎn)生的大規(guī)模的高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)中的2 類(lèi)敏感序列（短串聯(lián)重復(fù)序列和疾病相關(guān)序列）。在長(zhǎng)度為100～400 kbp 的長(zhǎng)序列數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果表明，并行算法CGPU-F3SR 在整體上獲得了較高的識(shí)別準(zhǔn)確率、查準(zhǔn)率和較低假陽(yáng)性率的同時(shí)，顯著提升了識(shí)別吞吐量。

近些年來(lái)，學(xué)者研究發(fā)現(xiàn)基因組測(cè)序數(shù)據(jù)中還存在“均聚物”錯(cuò)誤?？紤]到基因組數(shù)據(jù)錯(cuò)誤類(lèi)型復(fù)雜、數(shù)據(jù)規(guī)模巨大、數(shù)據(jù)敏感的特點(diǎn)，在本文研究工作基礎(chǔ)上，下一步的研究方向?qū)ǎ?) 研究設(shè)計(jì)能準(zhǔn)確識(shí)別出“均聚物”錯(cuò)誤的敏感序列識(shí)別算法；2) 探索設(shè)計(jì)新的數(shù)據(jù)索引結(jié)構(gòu)或借鑒其他索引數(shù)據(jù)結(jié)構(gòu)（如簡(jiǎn)潔de Bruign圖結(jié)構(gòu)、FM-index 結(jié)構(gòu)等），研究設(shè)計(jì)在CPU/GPU 混合體系結(jié)構(gòu)集群上實(shí)現(xiàn)的并行算法，以進(jìn)一步加速求解超大規(guī)模超長(zhǎng)序列的基因組數(shù)據(jù)敏感序列識(shí)別問(wèn)題；3)研究如何將本文提出的敏感序列識(shí)別并行算法與安全比對(duì)算法進(jìn)行級(jí)聯(lián)，以實(shí)現(xiàn)大規(guī)模的高錯(cuò)誤率長(zhǎng)序列基因組數(shù)據(jù)的安全比對(duì)。