孫 輝,鐘 誠

(廣西大學(xué) 計(jì)算機(jī)與電子信息學(xué)院,南寧 530004)

(廣西高校并行分布式計(jì)算技術(shù)重點(diǎn)實(shí)驗(yàn)室,南寧 530004)

1 引 言

人類基因組測序已經(jīng)在諸多領(lǐng)域取得了應(yīng)用.然而,如果基因組數(shù)據(jù)沒有得到有效的保護(hù),那么將會帶來很大的隱私泄露風(fēng)險[1-3].Gymrek等人的研究表明,通過短串聯(lián)重復(fù)序列(STRs)和單核苷酸多態(tài)性(SNPs),可以根據(jù)基因組數(shù)據(jù)和個人公開信息發(fā)起重識別攻擊(RIA),從而導(dǎo)致個人隱私數(shù)據(jù)泄露[4].基因組數(shù)據(jù)敏感序列識別是解決DNA數(shù)據(jù)隱私敏感問題的一個重要手段[5,6].

基因組敏感序列是人類DNA序列中易受隱私攻擊影響的核苷酸序列,分為短串聯(lián)重復(fù)序列和疾病相關(guān)序列兩類[6,7].短串聯(lián)重復(fù)序列是指基本重復(fù)單元為1～6bp的串聯(lián)重復(fù)序列[8].由于短串聯(lián)重復(fù)序列在人類遺傳變異中占比較大,且容易通過聚合酶鏈?zhǔn)椒磻?yīng)(PCR)擴(kuò)增,所以短串聯(lián)重復(fù)序列在疾病診斷和人員身份鑒定中具有重要作用.疾病相關(guān)序列是指攜帶疾病易感基因的核苷酸序列[7].由于疾病易感基因翻譯成蛋白質(zhì)可使個體表現(xiàn)出不同的性狀,所以疾病相關(guān)序列可被用于個體定位.所謂敏感序列識別是指從基因組原始測序數(shù)據(jù)中識別出以上兩類易受隱私攻擊影響的DNA敏感序列.通過敏感序列識別,對敏感序列片段給予保護(hù),可以避免個人隱私數(shù)據(jù)的泄露.

目前已有用于識別短串聯(lián)重復(fù)序列的一些啟發(fā)式算法REscan[9]、lobSTR[10]和mTR[11]等.然而,這些算法只能識別出部分基因組數(shù)據(jù)敏感序列,且對高錯誤率短串聯(lián)重復(fù)序列識別準(zhǔn)確性較差.為有效識別基因組數(shù)據(jù)敏感序列,Cogo等人提出了短序列過濾算法SRF[6].SRF算法的實(shí)現(xiàn)依賴于第三方數(shù)據(jù)庫收集的敏感序列和布隆過濾器(Bloom filter)[12,13].針對SRF算法拓展性差、難以識別潛在相似性敏感序列問題,Decouchant等人提出基于k-mer分割的長序列過濾算法LRF,并通過等位基因(allele)變異組合識別具有潛在相似性的敏感序列[7,14].Fernandes等人在LRF算法的基礎(chǔ)上,通過連鎖不平衡(LD)實(shí)現(xiàn)基因組數(shù)據(jù)敏感序列的多標(biāo)簽分類[15].Lambert等人通過敏感序列識別算法和軟件防護(hù)擴(kuò)展技術(shù)(SGX)實(shí)現(xiàn)了DNA序列的安全比對[16].隨著單分子實(shí)時測序技術(shù)(SMRT)和牛津納米孔測序技術(shù)(ONT)等第3代測序技術(shù)的發(fā)展,測序平臺產(chǎn)生的序列越來越長且具有高錯誤率的特點(diǎn)[17].

基于多哈希映射的布隆過濾器無法對潛在相似序列產(chǎn)生相同的哈希映射[12].若采用基于布隆過濾器的敏感序列識別算法對高錯誤率的基因組長序列進(jìn)行敏感序列識別,則其獲得的識別準(zhǔn)確率較低.為更有效地識別出基因組數(shù)據(jù)高錯誤率長序列中的敏感序列,本文提出一種融合過濾和相似度計(jì)算的敏感序列識別算法.本文剩余內(nèi)容組織如下:第2節(jié)將詳細(xì)闡述如何使用雙布隆過濾器避免對重復(fù)序列的多次計(jì)算以及如何改進(jìn)相似性度量計(jì)算方法以識別敏感序列;第3節(jié)給出本文算法和同類算法的實(shí)驗(yàn)評估結(jié)果;第4節(jié)總結(jié)本文工作.

2 方 法

設(shè)N[0:n-1]表示長度為n的第3代測序序列,本文融合過濾和相似度計(jì)算的敏感序列識別方法的過程是:第1步,分割序列N為n-r+1條短序列Rj=N[j:j+r-1],r為Rj的長度,j=0,1,2,…,n-r,n>r;第2步,采取雙布隆過濾器方法對短序列Rj進(jìn)行動態(tài)去重,以避免對相同短序列的重復(fù)計(jì)算,其中布隆過濾器BF1用于敏感序列去重,布隆過濾器BF2用于非敏感序列去重,j=0,1,…,n-r;第3步,對短序列Rj進(jìn)行短串聯(lián)重復(fù)序列識別,以判斷序列Rj是否是短串聯(lián)重復(fù)序列,若Rj是短串聯(lián)重復(fù)序列,則將Rj映射定位存儲至布隆過濾器BF1位數(shù)組B1,j=0,1,2,…,n-r;第4步,將短序列Rj與GWAS Catalog數(shù)據(jù)庫疾病敏感序列進(jìn)行計(jì)算比對,以判斷Rj是否是疾病相關(guān)序列,若Rj是疾病相關(guān)序列,則將Rj映射定位存儲至布BF1位數(shù)組B1,否則將Rj映射定位存儲至BF2位數(shù)組B2,j=0,1,2,…,n-r;第5步,根據(jù)短序列R0,R1,…,Rn-r的敏感識別結(jié)果res[0:n-1],生成序列N[0:n-1]的兩條掩碼序列Npub[0:n-1]和Npri[0:n-1]作為算法識別結(jié)果.圖1描述了本文方法的處理流程.

圖1 本文識別方法的處理流程

在圖1中,結(jié)果向量res[0:n-1]用于記錄每條短序列Rj的敏感序列識別中間結(jié)果;為保存算法最終識別結(jié)果,本文采用類似于文獻(xiàn)[7]的方法,根據(jù)結(jié)果向量res[0:n-1]生成序列N[0:n-1]的兩條掩碼序列Npub[0:n-1]和Npri[0:n-1],其中Npub表示掩碼非敏感序列,Npub中的敏感序列片段將以“*”掩碼,Npri表示掩碼敏感序列,Npri中的非敏感序列片段將以“*”掩碼.通過生成掩碼序列,對掩碼敏感序列進(jìn)行隱私保護(hù),以避免針對基因組數(shù)據(jù)的隱私攻擊.

下面詳細(xì)闡述融合過濾和相似度計(jì)算的高錯誤率敏感核苷酸序列識別方法.

2.1 序列分割

為識別測序長序列中的敏感序列,采取滑動窗口策略將長序列N[0:n-1]分割為n-r+1條、每條長度為r的短序列Rj=N[j:j+r-1],其中Rj和Rj+1共享重疊的r-1個核苷酸,以確保每條短序列敏感識別的準(zhǔn)確性,j=0,1,2,…,n-r,n>r.

人類DNA序列中,約有99.5%的序列具有高度一致性,只有0.5%的核苷酸記錄個體獨(dú)有信息[6].將長序列分割為短序列,會產(chǎn)生大量相同短序列.若對所有短序列均進(jìn)行敏感序列識別,則將耗費(fèi)許多計(jì)算資源.為避免對相同短序列的重復(fù)計(jì)算,本文利用布隆過濾器對分割得到的序列去重.

2.2 序列去重

布隆過濾器由一個二進(jìn)制向量和若干哈希函數(shù)組成[14].基于多哈希映射的布隆過濾器可以降低哈希沖突的影響.

2.2.1 初始化布隆過濾器

設(shè)布隆過濾器BF1和BF2分別用于敏感序列和非敏感序列去重,其位向量分別為B1[0:b1-1]和B2[0:b2-1],使用哈希函數(shù)的個數(shù)分別為h1和h2.初始化布隆過濾器,需計(jì)算位向量B1和B2的規(guī)模b1和b2[12]:

(1)

式(1)中,ti表示布隆過濾器BFi可容納序列數(shù)量的上限,pi表示誤報(bào)率,BFi位向量規(guī)模bi和可容納序列數(shù)量的上限ti呈正相關(guān),bi和誤報(bào)率pi(0

(2)

在計(jì)算出b1和b2,h1和h2之后,置B1,i=0、B2,k=0,i=0,1,2,…,b1-1,k=0,1,2,…,b2-1.

2.2.2 短序列去重

對短序列Rj去重,需要根據(jù)已知短序列Rk的識別結(jié)果將Rk映射定位存儲至對應(yīng)布隆過濾器,k=0,1,2,…,j-1.若Rk為敏感短序列,則對Rk進(jìn)行h1次哈希計(jì)算得到h1個哈希值H1,i(Rk),并將布隆過濾器BF1中位向量B1[H1,i(Rk)]置為1,i=0,1,2,…,h1-1.若Rk為非敏感短序列,則對Rk進(jìn)行h2次哈希計(jì)算得到h2個哈希值H2,i,并將布隆過濾器BF2中位向量B2[H2,i(Rk)]置為1,i=0,1,2,…,h2-1.

若要判斷Rj是否是一條重復(fù)序列,首先對Rj進(jìn)行h1次哈希計(jì)算,得到h1個哈希值H1,i(Rj),并在BF1中的位向量B1進(jìn)行查詢,若B1[H1,i(Rj)]均為1,i=0,1,2,…,h1-1,則Rj是一條重復(fù)短序列,且是敏感短序列,此時置res[j]=1.若B1[H1,i(Rj)]不全為1,則對Rj進(jìn)行h2次哈希計(jì)算,得到h2個哈希值H2,i(Rj),并在BF2中的位向量B2進(jìn)行查詢,若B2[H2,i(Rj)]均為1,i=0,1,2,…,h2-1,則Rj是重復(fù)短序列,且是非敏感短序列,此時置res[j]=0.若B2[H2,i(Rj)]不全為1,則說明Rj不存在于BF1或BF2中,Rj不是一條重復(fù)序列,需要對Rj執(zhí)行短串聯(lián)重復(fù)和疾病相關(guān)序列識別處理.

為更好體現(xiàn)算法去重思想,圖2給出序列去重示例.

圖2 布隆過濾序列去重示例

在圖2中,BF1和BF2的位向量B1和B2均是一個8位的二進(jìn)制向量;BFi序列去重使用3個哈希函數(shù)Hi,0、Hi,1和Hi,2,i=1,2;R1、R3和R4是敏感短序列,R2和R5是非敏感短序列,且R5是重復(fù)短序列(R5=R2),為此去重序列R5.首先,初始化布隆過濾器BF1和BF2的位向量B1、B2均為0,此時BF1和BF2均為空,沒有記錄任何短序列信息.然后,將R1～R4映射定位存儲至BF1和BF2的位向量B1和B2;此時B1第1,3,4,6,7位被設(shè)置為1,B2第2,3,6位被設(shè)置為1.最后,序列R5通過BF1進(jìn)行敏感序列去重時所產(chǎn)生的3個哈希標(biāo)記不全為1.因此需要經(jīng)過BF2進(jìn)行非敏感序列去重;此時對R5進(jìn)行哈希映射所產(chǎn)生的3個哈希值在BF2中標(biāo)記均為1,因此R5是重復(fù)短序列,且R5由BF2過濾出,所以R5是非敏感短序列.通過雙布隆過濾器處理可知序列R5為非敏感短序列,因此,不需要執(zhí)行后續(xù)短串聯(lián)重復(fù)序列和疾病相關(guān)序列識別操作.

通過對短序列執(zhí)行雙布隆過濾器去重操作,可避免對短序列的重復(fù)計(jì)算.若短序列不存在于布隆過濾器BF1或BF2,則需要對短序列進(jìn)行短串聯(lián)重復(fù)序列識別.

2.3 短串聯(lián)重復(fù)序列識別

設(shè)Σ={A,C,G,T}表示組成DNA序列的符號集,Σ+表示Σ上的非空字符串集,對于序列u∈Σ+,一個理想的短串聯(lián)重復(fù)序列可表示為[u]q,其中q表示序列u的重復(fù)次數(shù),|u|表示序列u的長度且1bp≤|u|≤6bp,u為最小基本重復(fù)單元[8].短串聯(lián)重復(fù)序列的內(nèi)部重復(fù)結(jié)構(gòu)表明,通過合理分割短串聯(lián)重復(fù)序列可以得到記錄短串聯(lián)重復(fù)序列基本重復(fù)單元的相同子段[11].為此,借鑒文獻(xiàn)[11]的思想,本文采用以下步驟判斷短序列是否是含有“插入”、“刪除”、 “替換”錯誤的短串聯(lián)重復(fù)序列:

第1步.等長度分割短序列Rj為4個長度為w的短片段,Wj={Wj,1,Wj,2,Wj,3,Wj,4},j=0,1,2,…,n-r;

第3步.計(jì)算短序列Rj局部片段Wj,2和Wj,3的相似度sim(Wj,2,Wj,3),并根據(jù)短串聯(lián)重復(fù)序列相似度閾值str_sim判斷Rj是否為短串聯(lián)重復(fù)序列,若Rj為短串聯(lián)重復(fù)序列,則置結(jié)果向量res[j]=1,j=0,1,2,…,n-r.

sim(Wj,2,Wj,3)=

(3)

圖3給出本文方法識別短串聯(lián)重復(fù)序列的示例.

圖3 短串聯(lián)重復(fù)序列識別示例

從圖3中可知,通過上述步驟計(jì)算出第2、3個短片段的相似性達(dá)92.79%,大于給定的判別閾值str_sim=0.64,因此R是一條短串聯(lián)重復(fù)序列.

通過對Rj中的兩個局部片段統(tǒng)計(jì)k-mer詞頻,可以提取Rj的k-mer統(tǒng)計(jì)特征.若Rj具有較少的“插入”、“刪除”、“替換”錯誤,通過k-mer編碼Rj的局部片段,可以確保局部片段詞頻具有較高的一致性,從而可以計(jì)算得到較高的相似度值,進(jìn)而識別Rj是否是具有潛在相似性的高錯誤率的短串聯(lián)重復(fù)序列,j=0,1,2,…,n-r.

2.4 疾病相關(guān)序列識別

通過對序列結(jié)構(gòu)分析處理可以識別出DNA序列中包含的短串聯(lián)重復(fù)序列.除了短串聯(lián)重復(fù)序列,DNA序列中還可能包含有疾病相關(guān)序列.為識別短序列Rj是否是疾病相關(guān)序列,本文通過構(gòu)建敏感序列詞典,檢索與Rj堿基含量相似的第三方數(shù)據(jù)庫中的疾病相關(guān)序列,采用皮爾遜相關(guān)系數(shù)對Rj和檢索到的疾病相關(guān)序列進(jìn)行相似度計(jì)算比對,以識別出高錯誤率疾病相關(guān)序列,j=0,1,2,…,n-r.

2.4.1 構(gòu)建敏感序列詞典

構(gòu)建敏感序列詞典的目的是為了降低大量短序列與疾病相關(guān)序列相似度計(jì)算的代價.本文使用二維數(shù)組DICTI[m][4k+3]表示敏感序列詞典,其中m表示從第三方數(shù)據(jù)庫提取的疾病相關(guān)序列的數(shù)量,k為采用k-mer編碼核苷酸序列的k值.構(gòu)建敏感序列詞典的方法是:

2.4.2 堿基含量相似性檢索

構(gòu)建敏感核苷酸序列詞典之后,對Rj進(jìn)行堿基含量相似性檢索,以進(jìn)一步降低Rj和敏感序列詞典中的序列相似度計(jì)算的工作量.文獻(xiàn)[11]的研究表明,潛在相似性序列具有堿基含量的一致相似性.序列Rj中G、C、A堿基含量GCARj=[GRj,CRj,ARj],本文采用式(4)的條件,在DICTI中檢索與Rj具有堿基含量相似的序列Di:

1-|GCARj[e]-GCADi[e]|≥dis_sim

(4)

式(4)中,i=0,1,2,…,m-1,j=0,1,2,…,n-r,e=0,1,2,dis_sim表示疾病相關(guān)序列識別相似度閾值.設(shè)queryB為通過序列堿基含量相似性檢索得到短序列Rj和Di滿足式(4)條件的Di下標(biāo)集合,則queryB=[q0,q1,…qv,…,qt],其中qv為DICTI中堿基含量相似序列的下標(biāo),t≤m-1.

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

(5)

當(dāng)序列中含有“插入”、“刪除”、“替換”錯誤時,采用k-mer編碼短序列,其差異只體現(xiàn)在變異核苷酸相鄰k個堿基處,通過合理取值k,可有效提取高錯誤率疾病相關(guān)序列堿基統(tǒng)計(jì)特征.將短序列Rj和第三方數(shù)據(jù)庫疾病相關(guān)序列進(jìn)行相似度計(jì)算,通過設(shè)置相似度閾值從而可以識別出和第三方數(shù)據(jù)庫具有較高相似性的疾病相關(guān)序列.

2.5 識別算法

設(shè)待識別長序列為N[0:n-1]、滑動窗口策略分割短序列Rj的窗口長度為r、雙布隆過濾器誤報(bào)率分別為p1和p2、初始化雙布隆過濾器可插入最大序列數(shù)目分別為t1和t2、k-mer編碼取值為k、短串聯(lián)重復(fù)序列識別相似度閾值為str_sim、疾病相關(guān)序列識別相似度閾值為dis_sim、參考基因組為HG、GWAS Catlog數(shù)據(jù)庫中疾病相關(guān)序列數(shù)據(jù)集為diseaseData(包含m組數(shù)據(jù)).

算法1形式描述了融合過濾和相似度計(jì)算的高錯誤率敏感序列識別算法(簡記為F3SR算法).

算法1. F3SR

輸入:N[0:n-1],r,p1,p2,t1,t2,k,str_sim,dis_sim,HG,diseaseData

輸出:序列N識別結(jié)果的兩條掩碼核苷酸序列Npri,Npub

Begin

1.依據(jù)p1,p2,t1和t2初始化布隆過濾器BF1和BF2位數(shù)組B1和B2;

2.依據(jù)diseaseData,r,HG和k構(gòu)建敏感序列詞典dicti[m,4k+3];

3.res[0:n-1]←0;

6.forj=0ton-rdo

7.Rj←N[j:j+r];//提取第j條短序列Rj

8. 對序列Rj執(zhí)行h1次哈希計(jì)算得到h1個哈希值H1,0(Rj),H1,1(Rj),…,H1,h1-1(Rj);

9.ifBF1中B1[H1,0(Rj)],B1[H1,1(Rj)],…,B1[H1,h1-1(Rj)]的值全為1then

10.res[j]←1;

11. 對序列Rj執(zhí)行h2次哈希計(jì)算得到h2個哈希值H2,0(Rj),H2,1(Rj),…,H2,h2-1(Rj);

12.elseifBF2中B2[H2,0(Rj)],B2[H2,1(Rj)],…,B2[H2,h2-1(Rj)]的值全為1then

13.res[j]←0;

14.else//編碼短序列局部片段

19.c←1;

20.else

21.c←0;

22.endif

24.ifsim(Wj,2,Wj,3)≥str_simthen

25.res[j]←1;

26. 將B1[H1,0(Rj)],B1[H1,1(Rj)],…,B1[H1,h1-1(Rj)]的值置為1;

27.else

28. 依據(jù)短序列Rj的GAC堿基含量對dicti進(jìn)行相似性序列檢索以獲得queryB;

29.ifqueryB=?then

30.res[j]←0;

31. 將B2[H2,0(Rj)],B2[H2,1(Rj)],…,B2[H2,h2-1(Rj)]的值置為1;

32.else

34.fori=0 to |queryB|-1do

37.ifsim(Rj,Di)≥dis-simthen

38.res[j]←1;

39. 將B1[H1,0(Rj)],B1[H1,1(Rj)],…,B1[H1,h1-1(Rj)]的值置為1;

40.goto(6);

41.endif

42.endfor

43.res[j]←0;

44. 將B2[H2,0(Rj)],B2[H2,1(Rj)],…,B2[H2,h2-1(Rj)]的值置為1;

45.endif

46.endif

47.endif

48.endfor

49.保存布隆過濾器BF1和BF2為二進(jìn)制文件;

50.Npub←N[0:n-1];//初始化掩碼非敏感序列

51.Npri←N[0:n-1];//初始化掩碼敏感序列

52.beginPos← 0;

53.whilebeginPos

54.ifres[beginPos]≠1then

55.beginPos←beginPos+1;

56.else

57.fori=beginPostobeginPos+rdo

58.ifres[beginPos+r-i]=1andi≠beginPos+rthen

59.forj=beginPostobeginPos+r-ido

60.Npub[j]←‘*’;

61.endfor

62.beginPos←beginPos+r-i;

63. goto(53);

64.elseifi=beginPos+rthen

65.forj=beginPostobeginPos+rdo

66.Npub[j]←‘*’;

67.endfor

68.beginPos←beginPos+r;

69. goto(53);

70.endif

71.endfor

72.endif

73.endwhile

74.forj=0ton-1do

//根據(jù)掩碼序列Npub生成掩碼敏感序列Npri

75.ifNpub[j]≠‘*’then

76.Npri[j]←‘*’;

77.endif

78.endfor

End.

算法F3SR步1初始化布隆過濾器時間為O(1);步2構(gòu)建敏感序列詞典所需時間為O(m×r×4k);步3初始化結(jié)果向量的時間為O(n);步4～步5計(jì)算布隆過濾器哈希函數(shù)個數(shù)時間為O(1);步34～步42最壞情況下所需的時間為O(m×4k);步6～步48時間為O(max(n×h1,n×h2,n×r×4k,n×4k,n×m×4k));步49保存二進(jìn)制文件時間為O(1);步50～步51初始化掩碼序列的時間為O(n);步52～步73生成掩碼非敏感序列的時間為O(n×r2);步74～步78生成掩碼敏感序列的時間為O(n).由于k、r、m、h1和h2均為某個常數(shù),所以算法時間復(fù)雜度為O(n).

算法F3SR通過k-mer編碼對DNA序列進(jìn)行詞頻統(tǒng)計(jì),將待識別序列中每個核苷酸的鄰接信息進(jìn)行有效提取.由于相似序列的差異只體現(xiàn)在差異核苷酸的局部序列片段,通過k-mer編碼可以使具有相似性的核苷酸序列保持詞頻的整體一致性,通過皮爾遜相關(guān)系數(shù)可計(jì)算出較高相似度得分的序列,所以可以識別出高錯誤率的敏感核苷酸序列.此外,基于“序列過濾”的思想,算法始終避免對重復(fù)序列的計(jì)算,有效提升了算法的效率.

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

3.1 實(shí)驗(yàn)環(huán)境與數(shù)據(jù)

實(shí)驗(yàn)在設(shè)立于廣西大學(xué)的國家高性能計(jì)算中心南寧分中心(2)https://hpc.gxu.edu.cn的Sugon 7000A超級并行計(jì)算機(jī)系統(tǒng)上進(jìn)行,使用的計(jì)算節(jié)點(diǎn)配置為2×Intel Xeon Gold 6230處理器、512GB內(nèi)存以及8×900GB外部存儲空間,運(yùn)行操作系統(tǒng)CentOS7.4.采用C++語言編程實(shí)現(xiàn)算法,布隆過濾器實(shí)現(xiàn)參考自ArashPartow(3)https://github.com/ArashPartow/bloom實(shí)現(xiàn)的混合hash版本.

本文采用GWAS Catalog數(shù)據(jù)庫[18]全基因組關(guān)聯(lián)研究數(shù)據(jù)集gwas_catalog_v1.0.tsv構(gòu)建敏感序列詞典,該數(shù)據(jù)集包含HIV-1 replication、Opioid sensitivity等4680種疾病,共計(jì)216250組數(shù)據(jù).每組數(shù)據(jù)包含疾病名稱DISEASE、染色體編號CHR_ID、染色體位置CHR_POS等34個標(biāo)識信息.實(shí)驗(yàn)數(shù)據(jù)采用的參考基因組為NCBI(4)https://www.ncbi.nlm.nih.gov開放下載的人類基因組HG38.實(shí)驗(yàn)所用的短串聯(lián)重復(fù)序列來自TRDB數(shù)據(jù)庫(5)https://tandem.bu.edu,疾病易感基因來自GWAS Catalog數(shù)據(jù)庫(6)https://www.ebi.ac.uk/gwas.參照文獻(xiàn)[6]的數(shù)據(jù)預(yù)處理規(guī)則,預(yù)處理每條序列長度為50bp的3個基準(zhǔn)數(shù)據(jù)序列集,每個基準(zhǔn)數(shù)據(jù)集的敏感序列和非敏感序列的數(shù)量比例為1∶7,實(shí)驗(yàn)數(shù)據(jù)信息如表1所示.

表1 實(shí)驗(yàn)數(shù)據(jù)集信息

為獲得含有錯誤信息的序列數(shù)據(jù)集,采用表1數(shù)據(jù)集作為基準(zhǔn)序列,通過隨機(jī)執(zhí)行“插入”、“刪除”、“替換”操作產(chǎn)生10組錯誤率為2%～20%的DNA序列進(jìn)行實(shí)驗(yàn).實(shí)驗(yàn)將本文算法F3SR與同類算法SRF[6]和LRF[7]進(jìn)行測試比較識別的準(zhǔn)確率(accuracy,acc)[6]、查準(zhǔn)率(precision,pre)[19]和假陽性率(falsepositiverate,fpr)[7]:

(6)

(7)

(8)

其中,TP表示算法識別出真實(shí)敏感序列的總數(shù)、TN表示識別出真實(shí)非敏感序列的總數(shù)、FN表示將敏感序列誤識別為非敏感序列的總數(shù)、FP表示將非敏感序列誤識別為敏感序列的總數(shù).準(zhǔn)確率反映算法的整體識別結(jié)果,準(zhǔn)確率越高算法識別效果越好.查準(zhǔn)率反映真實(shí)敏感序列和算法識別出的敏感序列的比率,查準(zhǔn)率越高算法越好.假陽性率反映數(shù)據(jù)中非敏感序列被誤識別為敏感序列的比率,假陽性率越低越好.

3.2 算法參數(shù)選取

為使得采用布隆過濾器對序列去重時具有較低的假陽性率,本文采用文獻(xiàn)[14]布隆過濾器實(shí)驗(yàn)參數(shù)pi=0.0001對布隆過濾器初始化,并設(shè)置布隆過濾器容納序列數(shù)量上限ti=104,當(dāng)實(shí)際插入序列的數(shù)目超過ti時,動態(tài)擴(kuò)充布隆過濾器的容納上限為當(dāng)前上限的2倍,i=1,2.

F3SR算法判別短序列Rj是否為短串聯(lián)重復(fù)序列時,需要對Rj第2和第3個短片段進(jìn)行k-mer編碼.通過對TRDB數(shù)據(jù)庫短串聯(lián)重復(fù)序列檢索,TRDB數(shù)據(jù)庫中短串聯(lián)重復(fù)序列長度最短為25bp[20],為使得F3SR算法可以準(zhǔn)確識別短串聯(lián)重復(fù)序列,F3SR算法取w=12bp.此時,滑動窗口策略分割短序列Rj的長度r=4w=48bp.

設(shè)err表示測序引入的錯誤率,第三代測序錯誤率高達(dá)10～20%[11],考慮測序過程引入最大錯誤率err=20%,算法相似度閾值t應(yīng)滿足t=1-max(err)=0.8,因此短串聯(lián)重復(fù)序列識別相似度閾值str-sim≤t2=0.802=0.64,疾病相關(guān)序列識別相似度閾值dis-sim=t=0.80.

F3SR算法通過k-mer編碼提取序列統(tǒng)計(jì)特征.為選取合適的參數(shù)k和相似性度量方法,本文通過實(shí)驗(yàn)測試k-mer編碼k取值、數(shù)據(jù)集序列錯誤率err和皮爾遜相關(guān)系數(shù)(pearson correlation coefficient)、斯皮爾曼等級相關(guān)系數(shù)(spearman′s rank correlation coefficient)、曼哈頓距離(manhattan distance)、肯德爾一致性系數(shù)(kendall′s consistency coefficient)對F3SR算法識別性能的影響.實(shí)驗(yàn)結(jié)果如圖4所示.

圖4 All-Together數(shù)據(jù)集上錯誤率、k值和相似性度量方法對算法F3SR的準(zhǔn)確率、查準(zhǔn)率和假陽性率的影響

在圖4中,橫坐標(biāo)表示數(shù)據(jù)集序列的錯誤率err(%),err=0表示序列未攜帶“噪聲數(shù)據(jù)”,err=20表示當(dāng)前基準(zhǔn)數(shù)據(jù)集中每條序列含有“插入”、“刪除”和“替換”錯誤的比例約為20%.圖4的結(jié)果表明:當(dāng)k≥4時算法采用pearson相關(guān)系數(shù)整體上具有最高的識別準(zhǔn)確率和識別查準(zhǔn)率,最低的識別假陽性率,其次是spearman等級相關(guān)系數(shù)和kendall一致性系數(shù),算法采用manhattan距離度量相似度時,整體上識別準(zhǔn)確率和查準(zhǔn)率最低,且具有較高的識別假陽性率.當(dāng)k取值4～6時,本文算法識別準(zhǔn)確率、查準(zhǔn)率和假陽性率趨于平穩(wěn).為使F3SR算法在保持較低識別假陽性率的同時具有較高識別準(zhǔn)確率和查準(zhǔn)率,選取k=5進(jìn)行k-mer編碼.

3.3 算法比對實(shí)驗(yàn)

為評估算法在不同錯誤率err的序列數(shù)據(jù)集上運(yùn)行的識別性能,實(shí)驗(yàn)測試了算法F3SR和SRF[6]、LRF[7]在All-STR數(shù)據(jù)集上運(yùn)行的準(zhǔn)確率acc、查準(zhǔn)率pre和假陽性率fpr,實(shí)驗(yàn)結(jié)果如表2所示.

表2 算法F3SR、SRF和LRF在All-STR數(shù)據(jù)集上運(yùn)行的準(zhǔn)確率、查準(zhǔn)率和假陽性率

表2實(shí)驗(yàn)結(jié)果表明,F3SR算法在數(shù)據(jù)集All-STR上運(yùn)行的識別準(zhǔn)確率均高于算法LRF和SRF,且數(shù)據(jù)集序列的錯誤率越低算法的識別準(zhǔn)確率越高,即使錯誤率達(dá)到20%時仍具有89.925%的較高準(zhǔn)確率.這是因?yàn)殄e誤率越低的短串聯(lián)重復(fù)序列在分割成序列片段時,越能保證局部片段的詞頻統(tǒng)計(jì)一致性,通過對局部片段進(jìn)行相似性度量便可以得到較高的局部相似度值,從而使得算法識別序列中短串聯(lián)重復(fù)序列的準(zhǔn)確率較高.F3SR算法在All-STR數(shù)據(jù)集上運(yùn)行的識別查準(zhǔn)率均高于SRF算法;當(dāng)All-STR數(shù)據(jù)集序列的錯誤率為4%～6%時,F3SR算法的查準(zhǔn)率略低于LRF算法;當(dāng)序列的錯誤率為2%、8%～20%時F3SR算法的查準(zhǔn)率高于LRF算法.這表明在高錯誤率的序列數(shù)據(jù)集上識別短串聯(lián)重復(fù)序列時,F3SR算法比其他兩個算法更有效.與算法LRF和SRF相比,F3SR算法在All-STR數(shù)據(jù)集上逆行的假陽性率最高;這是因?yàn)閷Χ绦蛄芯植科芜M(jìn)行相似性度量識別短串聯(lián)重復(fù)序列時,若短序列的局部片段具有極高相似性但并不是短串聯(lián)重復(fù)序列時,算法產(chǎn)生了錯誤判別.

在All-STR數(shù)據(jù)集上的實(shí)驗(yàn)表明,F3SR算法可以有效識別高錯誤率的短串聯(lián)重復(fù)序列.

為評估F3SR算法是否可以有效識別疾病相關(guān)序列,實(shí)驗(yàn)測試了算法F3SR、SRF[6]和LRF[7]在不同錯誤率err的序列數(shù)據(jù)集All-Disease上的識別性能,實(shí)驗(yàn)結(jié)果如表3所示.

表3 算法F3SR、SRF和LRF在All-Disease數(shù)據(jù)集上運(yùn)行的準(zhǔn)確率、查準(zhǔn)率和假陽性率

從表3可以看到,當(dāng)All-Disease數(shù)據(jù)集中序列錯誤率為2%～20%時:無論是某種具體錯誤率情形還是平均情況,F3SR算法整體上具有最高的識別準(zhǔn)確率和查準(zhǔn)率,LRF算法具有第2高的識別準(zhǔn)確率和查準(zhǔn)率,SRF算法的識別準(zhǔn)確率和查準(zhǔn)率最低;F3SR算法在All-Disease數(shù)據(jù)集上識別假陽性率最低,其次是SRF和LRF算法.

在All-Disease數(shù)據(jù)集上的實(shí)驗(yàn)表明,相較于SRF和LRF算法,F3SR算法識別疾病相關(guān)序列時,具有更高的準(zhǔn)確率和查準(zhǔn)率,具有最小的識別假陽性率,整體上識別效果較好.

為評估在既有短串聯(lián)重復(fù)序列又有疾病相關(guān)序列的數(shù)據(jù)集上運(yùn)行算法的識別性能,本文進(jìn)一步測試了3種算法在序列含有不同錯誤率err的真實(shí)混合數(shù)據(jù)集All-Together上的算法識別性能,實(shí)驗(yàn)結(jié)果如表4所示.

表4 算法F3SR、SRF和LRF在All-Together數(shù)據(jù)集上運(yùn)行的準(zhǔn)確率、查準(zhǔn)率和假陽性率

表4結(jié)果表明,對于各種錯誤率的序列,與算法SRF和LRF相比,F3SR算法在All-Together數(shù)據(jù)集上運(yùn)行時,整體上獲得最高的識別準(zhǔn)確率和查準(zhǔn)率.在假陽性率方面,對于錯誤率2%～12%的序列數(shù)據(jù)集,本文算法略低于LRF和SRF算法;對于錯誤率14%～20%的序列數(shù)據(jù)集,本文算法具有最低的假陽性率.算法F3SR相較于LRF和SRF,在錯誤率2%～20%的數(shù)據(jù)集All-Together平均識別準(zhǔn)確率分別提高1.96%和3.66%,平均查準(zhǔn)率分別提高40.08%和68.36%,平均假陽性率分別高了0.114%和0.038%.

綜上所述,相較于SRF和LRF算法,本文算法F3SR在3個數(shù)據(jù)集上運(yùn)行時,整體上具有更高的識別準(zhǔn)確率和識別查準(zhǔn)率,可有效識別具有潛在相似性的高錯誤率基因組數(shù)據(jù)敏感序列.

4 總 結(jié)

通過將長序列分割為多條短序列,對每條短序列采取雙布隆過濾方法進(jìn)行序列去重,可避免對每條短序列的重復(fù)計(jì)算;采用k-mer編碼策略對短序列進(jìn)行詞頻統(tǒng)計(jì)可提取具有潛在相似性序列的堿基統(tǒng)計(jì)特征,通過改進(jìn)序列相似度計(jì)算比對方法,可有效識別具有高錯誤率的敏感序列,使得識別算法獲得更高的識別準(zhǔn)確率和查準(zhǔn)率.融合長序列分割、過濾和相似度計(jì)算的方法可以有效識別出高錯誤率敏感序列,但算法需要花費(fèi)一些時間與第三方數(shù)據(jù)庫敏感序列進(jìn)行相似度計(jì)算比對,這對于序列長度越來越長的大規(guī)模測序數(shù)據(jù)集上運(yùn)行的時間開銷較高.如何設(shè)計(jì)并行計(jì)算方法加速長序列中敏感序列的識別過程,以適應(yīng)大規(guī)?；蚪M測序數(shù)據(jù)應(yīng)用,將是下一步的研究工作.