羅賢橦，鐘 誠，黎 瑤

(廣西大學(xué) 計(jì)算機(jī)與電子信息學(xué)院 廣西高校并行分布式計(jì)算技術(shù)重點(diǎn)實(shí)驗(yàn)室，南寧 530004)

1 引 言

以單分子實(shí)時(shí)測序技術(shù)(SMRT)[1]和Oxford Nanopore技術(shù)(ONT)[2]為代表的第三代測序平臺(tái)(Third Generation Sequencing，TGS)極大地降低了測序的成本[3]，可以解決人類基因組中具有復(fù)雜區(qū)域的測序問題[4].另一方面，第三代測序平臺(tái)產(chǎn)生錯(cuò)誤率超過15%、長度超過10k bp的長序列.測序序列長度的增加，使得一個(gè)長序列中可能包含一個(gè)完整的結(jié)構(gòu)變異區(qū)域，序列比對(duì)問題從僅需要處理較短長度的錯(cuò)誤(例如SNP錯(cuò)誤和較短的“indel”錯(cuò)誤)，演變到需要處理較長的結(jié)構(gòu)變異(轉(zhuǎn)置、易位、重復(fù)和長度超過50bp的“indel”)錯(cuò)誤，從而增加了序列比對(duì)的難度，降低了序列比對(duì)的敏感度.

針對(duì)這類錯(cuò)誤率較高且包含錯(cuò)誤類型較復(fù)雜的長序列比對(duì)方法有BWA-MEM[5]、BLASR[6]、Kart[7]、YAHA[8]、LAMSA[9]、Minimap2[10]和LordFAST[11]等.以BWA-MEM[5]、Kart[7]和Minimap2[10]為代表的長序列比對(duì)算法使用了雙端映射(paired-end mapping，PEM)的方法，它依賴于精確匹配將種子定位到結(jié)構(gòu)變異錯(cuò)誤的任意一側(cè)，且將種子擴(kuò)展到序列全部堿基.第三代測序數(shù)據(jù)中的高錯(cuò)誤率會(huì)導(dǎo)致比對(duì)算法難以找到長度較長的精確匹配，且序列長度的增長將會(huì)導(dǎo)致產(chǎn)生許多假陽性的較短的精確匹配結(jié)果，這將影響種子搜索的效率以及長序列比對(duì)的敏感度和查全率.隨著第三代測序技術(shù)的發(fā)展和序列組裝算法的改進(jìn)[12]，分割序列映射(split-read mapping，SRM)[8]成為長序列比對(duì)的一種有效策略.它更易于發(fā)現(xiàn)序列中的結(jié)構(gòu)變異錯(cuò)誤，且可在比對(duì)過程中減少結(jié)構(gòu)變異錯(cuò)誤對(duì)比對(duì)結(jié)果的影響.比對(duì)方法STAR[13]、YAHA[8]、LAMSA[9]、LordFAST[11]采用了該種策略.其思想是：首先將長序列分割成等長的若干片段，并將片段映射到參考基因組，生成片段定位的侯選位置，然后以片段侯選位置為依據(jù)，根據(jù)片段間可能存在的“錯(cuò)誤”關(guān)系將片段侯選位置依次連接起來，從而獲得長序列的比對(duì)結(jié)果.其中片段定位的過程將片段與參考基因比對(duì)，生成一組滿足給定編輯距離的片段侯選位置.這些比對(duì)算法在處理高錯(cuò)誤率的長序列比對(duì)時(shí)，查全率有所降低，且對(duì)錯(cuò)誤率高的真實(shí)數(shù)據(jù)集比對(duì)的敏感度較低.

為解決現(xiàn)有序列比對(duì)算法對(duì)長度更長、錯(cuò)誤率更高的第三代測序數(shù)據(jù)比對(duì)敏感度不夠高的問題，本文研究使用對(duì)高編輯距離閾值更敏感的基于Hash索引的變長種子擴(kuò)展方法，通過在種子驗(yàn)證階段特殊處理新型錯(cuò)誤“均聚物”[14]，以提高算法對(duì)第三代測序長序列比對(duì)的查全率，在連接序列片段前，加入對(duì)片段侯選位置的質(zhì)量評(píng)估，以過濾質(zhì)量不高的侯選位置、提高連接效率，且在連接片段候選位置時(shí)，對(duì)不同錯(cuò)誤類型賦予不同罰分，計(jì)算比對(duì)相似度，以減少假陽性的結(jié)果，獲得高準(zhǔn)確度和提高敏感度.

2 方 法

將高錯(cuò)誤率長序列與參考基因組進(jìn)行比對(duì)的第一步是將每條長序列劃分為若干較短的片段.第二步以基于Hash索引的變長種子擴(kuò)展算法定位序列片段在參考基因組中的候選位置，在種子驗(yàn)證過程中將連續(xù)“插入刪除”相同堿基錯(cuò)誤所帶來的編輯距離設(shè)置為1，使種子擴(kuò)展算法可以有效處理第三代測序長序列中的新類型錯(cuò)誤“均聚物”[14]，以使得片段定位時(shí)可提升查全率；采用全映射比對(duì)思想，尋找所有滿足編輯距離閾值的片段侯選位置.第三步連接序列片段的侯選位置，對(duì)片段侯選位置進(jìn)行評(píng)分，以去除質(zhì)量不高的侯選位置，然后將片段連接以構(gòu)建長序列比對(duì)結(jié)果的基本骨架.第四步，補(bǔ)全片段間的空隙，在補(bǔ)全過程中引入Z-drop[5]方法，以避免連接兩條不相關(guān)的片段而帶來的假陽性結(jié)果.本文給出的基于分割-全映射-篩選-連接-補(bǔ)全的高錯(cuò)誤率長序列比對(duì)方法的執(zhí)行過程如圖1所示.

下面詳細(xì)地闡述基于分割-全映射-篩選-連接-補(bǔ)全的高錯(cuò)誤率長序列比對(duì)方法.

2.1 片段定位

已有的種子擴(kuò)展方法被設(shè)計(jì)用于長度較短且錯(cuò)誤率較低的短序列比對(duì)(short-read alignment)[15].這類種子擴(kuò)展方法一般采用BWT索引結(jié)構(gòu)[16]或FM索引結(jié)構(gòu)[17]儲(chǔ)存參考基因組序列信息，這兩種索引結(jié)構(gòu)沒有存儲(chǔ)全部的序列信息.因此，對(duì)于高編輯距離閾值的比對(duì)敏感度較低.為了提升對(duì)高編輯距離閾值的序列比對(duì)的敏感度，本文采用對(duì)編輯距離閾值更加敏感的Hash索引結(jié)構(gòu)[18]，通過儲(chǔ)存參考基因組的k階子串k-mers出現(xiàn)的所有位置，來實(shí)現(xiàn)快速搜索大量的序列[22].與BWT[16]、FM-index[17]等索引機(jī)制相比，Hash索引對(duì)高編輯距離閾值和長度更長的序列比對(duì)具有更高的敏感度，可以找到更多的符合閾值的比對(duì)結(jié)果.但隨著編輯距離閾值的增大，會(huì)出現(xiàn)更多的假陽性的侯選位置，此外，Hash索引結(jié)構(gòu)的不足是難以處理序列中存在的“插入刪除(indel)”錯(cuò)誤.為解決上述問題，本文采用變長種子選擇算法[19,20]處理“插入刪除(indel)”錯(cuò)誤，且在種子驗(yàn)證過程中，將連續(xù)“插入刪除”多個(gè)相同堿基的錯(cuò)誤的編輯距離設(shè)置為1，以使得在種子擴(kuò)展時(shí)可以處理第三代測序長序列存在的新型“均聚物”錯(cuò)誤[14]，同時(shí)采用附加k-mer算法[21]去除可能產(chǎn)生的假陽性錯(cuò)誤.

圖1 基于分割-全映射-篩選-連接-補(bǔ)全的高錯(cuò)誤率長序列比對(duì)過程Fig.1 Process of aligning long reads with high error rate based on split-all mapping-filter-linking-completion

2.1.1 變長種子生成

與第二代測序數(shù)據(jù)主要存在“替換”錯(cuò)誤相比，第三代測序數(shù)據(jù)中存在的錯(cuò)誤類型還有“插入刪除”錯(cuò)誤，以往使用漢明距離確定種子候選位置的方式，并不適用于第三代測序數(shù)據(jù).為使得比對(duì)算法既支持處理“替換”錯(cuò)誤，又支持處理“插入刪除”錯(cuò)誤，本文采用可以處理“插入刪除”錯(cuò)誤的變長種子篩選方法[19,20]，并執(zhí)行附加k-mer算法[21]以去除假陽性的侯選位置.

設(shè)測序序列片段和參考基因組之間允許的編輯距離為e，根據(jù)鴿巢原理[5]，則至少選取序列片段的e+1個(gè)非重疊的k-mers，才能確保至少有一個(gè)k-mer與參考基因組的精確匹配是正確的.因此，通過將一組e+1個(gè)不重疊的固定長度的k-mer與參考基因組比對(duì)，即可獲得一組侯選位置集合Seed.為了使得算法可以處理“插入刪除”錯(cuò)誤，采用充分利用不重疊的k-mer之間堿基的變長種子算法，通過將固定長度k-mer進(jìn)行擴(kuò)展，獲得一組新的變長種子侯選位置集合BSeed.在擴(kuò)展時(shí)為了提高侯選位置驗(yàn)證階段的效率，優(yōu)先擴(kuò)展出現(xiàn)頻率較高的k-mer，以降低變長種子候選集BSeed中候選位置的數(shù)量.經(jīng)過這樣處理，可以最大程度地利用待比對(duì)片段上的所有堿基進(jìn)行種子的生成，使得算法可以處理“插入刪除”錯(cuò)誤，提高片段比對(duì)的準(zhǔn)確度.圖2展示了一組3-mer的變長種子的生成過程實(shí)例，其中圖2(a)展示找到一組不重疊的固定長度的k-mer與參考基因組進(jìn)行精確匹配所獲得的一組侯選位置，其中3-mer GGG在參考基因組中有5個(gè)侯選位置.圖2(b)展示對(duì)固定長度的k-mer進(jìn)行擴(kuò)展，通過將3-mer GGG擴(kuò)展為GAGGG，得到長度為5的精確匹配，獲得一組假陽性較少的侯選位置.

圖2 變長種子生成過程實(shí)例Fig.2 Example of generating variable-length seeds

同時(shí)，為增加種子定位的準(zhǔn)確度，采用附加k-mer的方法[21]，將選取比對(duì)片段的e+2個(gè)非重疊的k-mers，作為一組候選種子，以確保其中至少有2個(gè)k-mer與參考基因組的精確匹配是正確的，進(jìn)一步減少假陽性的種子侯選位置.

2.1.2 種子驗(yàn)證

在確定候選位置之后，對(duì)種子候選位置進(jìn)行驗(yàn)證，僅保留所有滿足編輯距離閾值小于等于e的侯選位置，作為最終的片段位置.第三代測序平臺(tái)產(chǎn)生的測序數(shù)據(jù)中可能含有新型的 “均聚物”錯(cuò)誤[16]，即出現(xiàn)連續(xù)“插入或刪除多個(gè)相同堿基”的錯(cuò)誤.當(dāng)遇到這種錯(cuò)誤時(shí)，驗(yàn)證操作會(huì)將原本定位正確的種子過濾掉，從而降低比對(duì)的查全率.針對(duì)這種新類型的錯(cuò)誤，本文方法在驗(yàn)證種子時(shí)，將“連續(xù)插入或刪除同一個(gè)堿基的錯(cuò)誤”作為一個(gè)編輯錯(cuò)誤處理，即在驗(yàn)證種子階段，當(dāng)發(fā)現(xiàn)連續(xù)出現(xiàn)相同堿基的“插入/刪除“錯(cuò)誤時(shí)，則將其編輯距離定義為1，即將該種錯(cuò)誤所造成的不匹配作為1個(gè)編輯操作，以增加片段定位的查全率，從而提升長序列比對(duì)的敏感度.

2.2 骨架構(gòu)建

2.2.1 篩選片段侯選

在種子定位過程中，通過增大編輯距離來增加序列片段的比對(duì)侯選位置數(shù)量，這將導(dǎo)致假陽性侯選位置增多.因此，在骨架構(gòu)建之前，對(duì)片段侯選位置進(jìn)行篩選，以去除質(zhì)量不高的序列片段比對(duì)的結(jié)果.

對(duì)于含有許多重復(fù)的相同堿基的序列片段，將其定位到同樣具有大量重復(fù)相同堿基的參考基因組(例如人類基因組)位置上，會(huì)導(dǎo)致將片段定位到參考基因組中多個(gè)基因的多個(gè)不同位置上.這種片段定位得到的比對(duì)結(jié)果質(zhì)量不高，會(huì)造成最終比對(duì)的假陽性結(jié)果的上升，從而影響最終比對(duì)準(zhǔn)確度.通過統(tǒng)計(jì)長序列分割出的片段侯選位置，可以計(jì)算得到片段Pi的比對(duì)質(zhì)量：

(1)

其中sumi為片段Pi的全部比對(duì)結(jié)果數(shù)量，sij為片段Pi在參考基因組中第j條參考基因上的比對(duì)結(jié)果數(shù)量.當(dāng)Qi小于給定閾值時(shí)，舍棄該片段在第j條參考基因的比對(duì)結(jié)果.

2.2.2 片段連接

根據(jù)片段間的相互關(guān)系，將篩選后的片段進(jìn)行連接，以進(jìn)一步減少片段定位候選區(qū)域中的假陽性區(qū)域.從序列P0開始依次將片段的侯選位置Mi與其之前所有已連接的片段相連接,選取其中連接得分最高的連接結(jié)果，作為Mi的連接得分f(i)，直至所有侯選位置均被連接到至少一條骨架上為止.本文改進(jìn)后的連接得分f(i)的計(jì)算公式如下：

(2)

其中scorei為片段的侯選位置Mi與參考基因組的匹配得分，f(j)是片段Pj之前的j-1條片段的候選位置的連接得分，a(j,i)是片段Pi與它相連的片段Pj之間的關(guān)系罰分.測序序列共存在4種錯(cuò)誤情況(如圖3所示)：

4)當(dāng)Mi與Mj在參考基因組上的匹配方向相反時(shí)，片段Mi與Mj之間存在轉(zhuǎn)置錯(cuò)誤；其中ε是用于計(jì)算“插入刪除”錯(cuò)誤長度的參數(shù),τ,θ是用于計(jì)算給定“復(fù)制”結(jié)構(gòu)變異錯(cuò)誤長度的參數(shù)[9].

圖3 測序序列的4種錯(cuò)誤類型Fig.3 Four error types of sequencing sequences

將所有測序序列片段侯選位置連接生成一個(gè)有向無環(huán)圖(DAG)，從圖的起點(diǎn)到終點(diǎn)的一條路徑，即為一條長序列比對(duì)的連接結(jié)果，將生成的所有連接的得分進(jìn)行排序，找到連接得分較大排序在前面的那些連接結(jié)果.

2.3 空隙填補(bǔ)

當(dāng)片段連接完成后，對(duì)每條長序列的前幾個(gè)連接結(jié)果執(zhí)行最終比對(duì).為了完成候選區(qū)域的堿基到堿基的對(duì)準(zhǔn)，它以連接得分從高到低將連接好的片段之間的空隙進(jìn)行匹配，匹配過程使用動(dòng)態(tài)規(guī)劃全局比對(duì)方法[23]，將序列片段間的堿基一一對(duì)應(yīng)到參考基因組侯選位置上.設(shè)定閾值，當(dāng)兩個(gè)片段對(duì)準(zhǔn)得分(匹配為0，不匹配為-1)超過閾值時(shí)，則認(rèn)為該條連接結(jié)果無效，這樣減少將兩個(gè)不相關(guān)的片段連接到一起而得到假陽性的連接結(jié)果.

長序列比對(duì)算法主要由兩大部分工作組成，第一部分工作為分割序列及片段定位，通過將長序列分割為若干片段，以更快地得到較短的精確匹配的種子，對(duì)種子進(jìn)行驗(yàn)證，得到所有滿足編輯距離閾值的片段侯選位置.第二部分工作首先通過片段質(zhì)量評(píng)分，去除質(zhì)量不高的片段候選位置，通過片段間的連接關(guān)系去除假陽性片段侯選位置，以確保比對(duì)結(jié)果的準(zhǔn)確度.詳細(xì)比對(duì)過程如算法1所示.

算法1.基于分割-全映射-過濾-連接-補(bǔ)全的比對(duì)算法

輸入：編輯距離閾值e，測序序列ri,i=0,1,…,rnum

輸出：序列ri比對(duì)結(jié)果Si,i=,1,…,rnum

1.Begin

2.fori=1 tornumdo

3. 將序列ri以固定間距l(xiāng)en分割為n個(gè)片段P1,…,Pn;

4.forj=1tondo

//找到所有滿足編輯距離閾值e的片段候選位置選擇

//Pj的e+2個(gè)非重疊k-mer作為變長種子侯選集Seed；

5.Seed[1].start←1;//種子的起始位置

6.Seed[e].end←|Pj|-1;//種子的結(jié)束位置

7.fork=1toe+2do

//將e+2個(gè)種子中出現(xiàn)最多的種子擴(kuò)展，f為種子出現(xiàn)數(shù)

8.ifSeed[k].f>=Seed[k-1].fthen

9.Seed[k].start←Seed[k-1].end+1;

10.else

11.Seed[k-1].end←Seed[k].start-1;

12.endif

13.endfor

14.BSeed←Seed[1:e+2];

//將擴(kuò)展后的種子Seed[1..e+1]作為新侯選集BSeed

//在Hash索引中搜索BSeed中所有種子的位置，將其存

//入片段侯選位置List[1..num]中

15.forl=1tonumdo

//驗(yàn)證侯選位置，篩選出符合編輯距離e的侯選位置

16.edit←0；

17.forp=1tolendo//驗(yàn)證

18.if片段中第p位堿基與參考基因中List[l]+p和List[l]+p-1位堿基均不同then

19.edit←edit+1；

20.endif

21.endfor

22.endfor

23. Ifedit<=ethen//編輯距離小于閾值

24.scorei←|Pj|-edit;

26.endif

27.endfor

28.fork=1 tondo

//統(tǒng)計(jì)Pk在參考基因j上的候選位置數(shù)量sij

29.sumk←0;//篩選過程sumk為片段候選結(jié)果

30.forj=1 to mdo

31.sumk←sk+sij;

32.endfor

33.ifskj/sumk

34. 舍棄Pk中所有在參考基因j上的比對(duì)結(jié)果

35.endif

36.endfor

37.forl=1 tomdo

//遍歷所有片段，連接侯選位置Mi生成無環(huán)圖

//計(jì)算候選位置Ml的最佳連接得分

39.endfor

40.根據(jù)連接得分將ri的所有連接結(jié)果進(jìn)行排序；

41.選擇得分最高的前幾個(gè)連接結(jié)果進(jìn)行片段間的空隙填補(bǔ)，計(jì)算序 列ri與參考基因組的匹配得分；

42.將填補(bǔ)后的序列匹配得分排序；

43.將最高匹配得分的結(jié)果序列在參考基因組上的起始位置，作為序 列ri最終的比對(duì)結(jié)果集合S；

44.endfor

45.end

與同類的序列比對(duì)算法LAMSA[9]和LordFAST[11]相比,本文算法采用對(duì)高編輯距離更敏感的Hash索引結(jié)構(gòu)存儲(chǔ)參考基因組，在定位序列分割片段時(shí)，通過變長種子算法[13]使得基于Hash索引的種子擴(kuò)展可以處理“插入刪除(indel)”錯(cuò)誤的編輯距離，且將“連續(xù)插入/刪除多個(gè)相同堿基”錯(cuò)誤作為1個(gè)編輯距離，處理第三代測序序列特有的新型“均聚物”錯(cuò)誤，以獲得更全的片段侯選位置，確保片段定位的查全率，根據(jù)片段間的關(guān)系過濾掉片段定位侯選位置中的假陽性結(jié)果，在確保比對(duì)準(zhǔn)確度的同時(shí)，提升查全率和敏感度.

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

實(shí)驗(yàn)使用的計(jì)算機(jī)為4核Intel(r)Xeon CPU E5-2600 V2處理器、內(nèi)存容量128GB.運(yùn)行操作系統(tǒng)Ubuntu 16.04.采用C語言編程實(shí)現(xiàn)算法.通過網(wǎng)站NCBI(1)https://www.ncbi.nlm.nih.gov/home/download/下載了公共的人類基因組數(shù)據(jù)集hg19進(jìn)行比對(duì)實(shí)驗(yàn).

3.1 模擬數(shù)據(jù)實(shí)驗(yàn)

模擬實(shí)驗(yàn)數(shù)據(jù)采用Wgsim(2)https://github.com/lh3/wgsim以人類基因組hg19作為參考基因組生成錯(cuò)誤率分別為5%和10%的模擬序列數(shù)據(jù)，采用PBsim[23]以hg19作為參考基因組生成錯(cuò)誤率分別為15%、20%和25%的含有“插入刪除”和結(jié)構(gòu)變異錯(cuò)誤的長度大于等于5000bp的模擬序列.

對(duì)于生成的每條模擬序列，Wgsim和PBsim[23]都提供了模擬序列在參考基因組上的映射位置.一個(gè)模擬序列在參考基因組上的“真正”映射位置是已知的，當(dāng)比對(duì)序列得出的比對(duì)位置與其模擬映射位置差值在30bp以內(nèi)，則認(rèn)定為該序列被比對(duì)到參考基因組中正確的位置上[7].

序列比對(duì)準(zhǔn)確度precision和查全率recall的計(jì)算[24]：

(3)

(4)

其中TP為正確映射到參考基因上的序列數(shù)量,RN為被映射到參考基因組的序列數(shù)量，N為參加映射比對(duì)的序列數(shù)量.

在模擬數(shù)據(jù)上的實(shí)驗(yàn)主要是為了評(píng)估序列比對(duì)算法的準(zhǔn)確度和查全率.隨著測序序列錯(cuò)誤率的上升，比對(duì)算法采用的編輯距離也需隨之增加，以獲得更多的侯選位置.為此，首先對(duì)不同錯(cuò)誤率下編輯距離閾值的選擇進(jìn)行實(shí)驗(yàn)：對(duì)10000條長度為5000bp的錯(cuò)誤率分別為5%、10%、15%、20%和25%的模擬序列數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，測試不同的編輯距離閾值對(duì)算法的比對(duì)準(zhǔn)確度和查全率的影響.錯(cuò)誤率低于15%的模擬序列屬于低錯(cuò)誤率數(shù)據(jù)，選取編輯距離閾值小于5進(jìn)行實(shí)驗(yàn)，結(jié)果如表1所示.而錯(cuò)誤率大于等于15%的模擬序列屬于高錯(cuò)誤率數(shù)據(jù)，需要更高的編輯距離閾值，選取編輯距離閾值為5、10和15進(jìn)行對(duì)比試驗(yàn)，結(jié)果如表2所示.

表1和表2的實(shí)驗(yàn)結(jié)果表明，測序序列的錯(cuò)誤率越高，片段中所含的錯(cuò)誤個(gè)數(shù)也就越多.此時(shí)，如果不提高編輯距離的閾值，那么將會(huì)導(dǎo)致丟失許多的真陽性位置，從而導(dǎo)致比對(duì)準(zhǔn)確度和查全率的下降.當(dāng)測序序列中錯(cuò)誤率小于15%時(shí)，比對(duì)準(zhǔn)確度變化不大，考慮算法處理時(shí)間，編輯距離閾值選擇2即可.當(dāng)錯(cuò)誤率大于等于15%時(shí)，編輯距離閾值應(yīng)選擇15以獲得更高的比對(duì)準(zhǔn)確度和查全率.

表1 不同編輯距離閾值、低錯(cuò)誤率模擬長序列比對(duì)的準(zhǔn)確度與查全率Table 1 Accuracy and recall of simulated long-read alignment with different edit distance threshold and low error rate

表2 不同編輯距離閾值、高錯(cuò)誤率模擬長序列比對(duì)的準(zhǔn)確度與查全率Table 2 Accuracy and recall of simulated long-read alignment for different edit distance threshold and high error rate

下面，通過對(duì)模擬序列數(shù)據(jù)進(jìn)行兩組實(shí)驗(yàn)以評(píng)估不同序列比對(duì)算法的性能.

第一組比對(duì)實(shí)驗(yàn)：對(duì)長度相同、錯(cuò)誤率不同的模擬測序序列進(jìn)行實(shí)驗(yàn)，測試了LAMSA、LordFAST和本文算法HSSM的比對(duì)準(zhǔn)確度和查全率.LAMSA和LordFAST算法中閾值參數(shù)選擇其論文中的給定值(錯(cuò)誤率小于15%時(shí),容忍錯(cuò)誤率為4%，錯(cuò)誤率大于等于15%，容忍錯(cuò)誤率30%)，本文算法HSSM對(duì)于錯(cuò)誤率分別為5%和10%的序列比對(duì)編輯距離閾值取值為2，對(duì)于錯(cuò)誤率分別為15%、20%和25%的序列比對(duì)編輯距離閾值取值為15.實(shí)驗(yàn)結(jié)果如表3所示.

表3 不同錯(cuò)誤率下算法比對(duì)10000條長度為5000bp的序列準(zhǔn)確度和查全率Table 3 Accuracy and recall rate of aligning 10,000 longreads with length of 5000bp for different error rates

表3結(jié)果表明，算法HSSM采用了對(duì)高編輯距離閾值更加敏感的Hash索引方式，且在片段定位過程中針對(duì)第三代測序序列中的新錯(cuò)誤類型“均聚物”進(jìn)行處理，將所有滿足編輯距離閾值的片段比對(duì)結(jié)果作為片段連接時(shí)的侯選位置，整體上獲得更高的比對(duì)準(zhǔn)確度，且獲得了更多的片段候選位置，提升了定位片段的查全率，從而使得最終比對(duì)結(jié)果的查全率也隨之提高.

第2組比對(duì)實(shí)驗(yàn)選取錯(cuò)誤率分別為5%、10%、15%、20%和25%的10000條序列.對(duì)于錯(cuò)誤率為5%和10%的情形，選取長度分別為1000bp、2000bp、5000bp和10000bp的4組模擬序列進(jìn)行實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果見表4和表5.對(duì)于錯(cuò)誤率為15%、20%和25%的情形，由于是模擬第三代測序的長序列數(shù)據(jù)，所以選取長度為5000bp和10000bp的模擬長序列進(jìn)行對(duì)比，結(jié)果如表6、表7和表8所示.

從表4-表8的結(jié)果看，與LAMSA和LordFAST算法相比，在錯(cuò)誤率相同且錯(cuò)誤率較低的情形下，總體而言，本文算法HSSM的準(zhǔn)確度和查全率有所提升.隨著錯(cuò)誤率的升高，HSSM算法獲得的比對(duì)查全率高于其他兩種算法.這是由于本文算法HSSM采用了基于Hash索引結(jié)構(gòu)的變長種子算法，更加適應(yīng)高編輯距離閾值的序列比對(duì)，可以獲得更多的侯選位置，提升了片段定位的查全率，進(jìn)而提升了長序列比對(duì)的查全率，且通過侯選位置的篩選和片段間可能存在的錯(cuò)誤關(guān)系的罰分，有效去除部分假陽性錯(cuò)誤的結(jié)果，使得長序列比對(duì)保持高的準(zhǔn)確度.另一方面，隨著模擬序列長度的增長，將長序列分割成了更多的片段，獲得了更多的侯選位置，可以有效地根據(jù)片段侯選位置之間的位置關(guān)系，去除錯(cuò)誤的比對(duì)結(jié)果，使得算法準(zhǔn)確度隨著測序序列長度的上升而提高.

表4 錯(cuò)誤率為5%的10000條序列比對(duì)的準(zhǔn)確度和查全率Table 4 Accuracy and recall of aligning 10,000 long readswith error rate 5%

表5 錯(cuò)誤率為10%的10000條序列比對(duì)的準(zhǔn)確度和查全率Table 5 Accuracy and recall of aligning 10,000 long readswith error rate 10%

表6 錯(cuò)誤率為15%的10000條長序列比對(duì)準(zhǔn)確度和查全率Table 6 Accuracy and recall of aligning 10,000 long readswith error rate 15%

表7 錯(cuò)誤率為20%的10000條長序列比對(duì)準(zhǔn)確度和查全率Table 7 Accuracy and recall of aligning 10,000 long readswith error rate 20%

表8 錯(cuò)誤率為25%的10000條長序列比對(duì)準(zhǔn)確度和查全率Table 8 Accuracy and recall of aligning 10,000 long readswith error rate 25%

3.2 真實(shí)數(shù)據(jù)上的實(shí)驗(yàn)

對(duì)于真實(shí)數(shù)據(jù)集上的實(shí)驗(yàn)，采用比對(duì)敏感度sensitivity來評(píng)估序列比對(duì)算法[24]：

(5)

其中N為匹配對(duì)準(zhǔn)的序列條數(shù)，RN為測序數(shù)據(jù)集的序列總數(shù).

采用PacBio測序平臺(tái)產(chǎn)生的M130929數(shù)據(jù)集中的幾組真實(shí)數(shù)據(jù)進(jìn)行測試實(shí)驗(yàn)，結(jié)果如表9所示.

表9 真實(shí)數(shù)據(jù)集上算法的比對(duì)敏感度(%)Table 9 Sensitivity (%)of alignment algorithms on real dataset

表9的結(jié)果表明，算法HSSM使用了對(duì)高錯(cuò)誤率更加敏感的基于hash索引結(jié)構(gòu)的變長種子算法，使得hash索引可以處理“插入/刪除”的錯(cuò)誤類型，并針對(duì)第三代測序序列中新出現(xiàn)的“均聚物”錯(cuò)誤類型設(shè)定其編輯距離為1進(jìn)行比對(duì)，獲得了更多片段候選區(qū)域，從而獲得更多的比對(duì)結(jié)果，且通過采用附加k-mer算法減少了產(chǎn)生假陽性的比對(duì)結(jié)果，進(jìn)而獲得更高的比對(duì)敏感度.

綜合在模擬和真實(shí)序列數(shù)據(jù)上的實(shí)驗(yàn)結(jié)果表明，本文算法HSSM針對(duì)第三代測序數(shù)據(jù)的特點(diǎn)設(shè)計(jì)改進(jìn)，在獲得更多侯選位置的同時(shí)，去除了假陽性的比對(duì)結(jié)果，既確保獲得高的比對(duì)準(zhǔn)確度，又獲得了更高的查全率和敏感度.

4 總 結(jié)

通過基于Hash索引結(jié)構(gòu)的變長種子算法，采用全映射思想定位長序列分割出的序列片段，可以最大程度地確保分割出的序列片段比對(duì)的查全率.分割映射的序列比對(duì)方法可以更好地處理第三代測序數(shù)據(jù)中的結(jié)構(gòu)變異帶來的“插入/刪除”和“均聚物”錯(cuò)誤,且通過篩選和動(dòng)態(tài)連接片段侯選位置，可以得到更高的比對(duì)準(zhǔn)確率和敏感度.高錯(cuò)誤率的長序列比對(duì)算法中，編輯距離閾值的提高會(huì)增加種子候選數(shù)量，且隨著第三代測序序列的長度越來越長，也會(huì)導(dǎo)致種子候選位置增加，從而使得種子驗(yàn)證成本增加，算法的時(shí)間開銷也隨之上升.第三代測序數(shù)據(jù)的高錯(cuò)誤率的長序列比對(duì)是一個(gè)計(jì)算復(fù)雜問題.下一步工作將在借鑒全映射比對(duì)并行算法[25]的基礎(chǔ)上，研究設(shè)計(jì)求解高錯(cuò)誤率長序列比對(duì)問題的CPU/GPU混合并行算法，以確保比對(duì)準(zhǔn)確度、查全率和敏感度的同時(shí)，顯著加速比對(duì)完成.