增強(qiáng)型稀疏后綴數(shù)組索引的高錯(cuò)誤率reads比對(duì)

韋 好，鐘 誠(chéng)

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

E-mail:394178293@qq.com

1 引言

將高通量DNA測(cè)序儀產(chǎn)生的read序列映射(比對(duì))到參考基因組是生物信息學(xué)中十分重要的研究問(wèn)題之一.映射(Mapping)過(guò)程是高通量測(cè)序技術(shù)(HTS)數(shù)據(jù)處理中的第一步也是最耗時(shí)的步驟，它采用索引串匹配的方法將read序列與參考基因組進(jìn)行比對(duì)，確定每一條read在參考基因組上的位置，這一過(guò)程在基因組分析流程中起著關(guān)鍵作用.序列比對(duì)算法大致分為兩大類:基于哈希表索引的算法和基于前綴/后綴樹(shù)(數(shù)組)索引的算法［1］.第一類算法通常為參考基因組建立一個(gè)哈希索引表，根據(jù)“種子(seed)-擴(kuò)展”策略，將比read短得多的子序列(稱為種子)直接定位到參考基因組，然后根據(jù)每個(gè)映射位置上的種子擴(kuò)展的結(jié)果來(lái)確定read的位置.這類算法雖能有效檢測(cè)識(shí)別出單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)或結(jié)構(gòu)變異，但其最大的缺陷是內(nèi)存空間占用很大.基于哈希表索引的代表性的比對(duì)算法有BFAST［2］、SOAP［3］、MAQ［4］、YAHA［5］和 FSVA［6］.第二類算法通常搜索參考基因組的前綴/后綴樹(shù)(數(shù)組)，然后利用BWT(Burrows-Wheeler Transform)［7］方法計(jì)算每條 read 在參考基因組中的位置.這類比對(duì)算法可以節(jié)省大量存儲(chǔ)空間和有效地減少參考序列與read序列之間的非精確匹配，但其比對(duì)速度較慢.基于前綴/后綴樹(shù)(數(shù)組)的代表性比對(duì)算法有BWA［8］、Bowtie2［9］、SOAP2［10］和 HPG［11］.

在單分子測(cè)序儀的最新進(jìn)展中，研究人員獲得了比Sanger測(cè)序儀提供的長(zhǎng)度更長(zhǎng)的reads.自2010年，Pacific Biosciences發(fā)布的第一個(gè)實(shí)時(shí)單分子測(cè)序儀PacBio RS以來(lái)，單分子測(cè)序儀產(chǎn)生的reads長(zhǎng)度一直在增加［12］.另一方面，long read(長(zhǎng)序列)包含參考基因組中較高的錯(cuò)誤率(約為15%)、更多結(jié)構(gòu)變異以及錯(cuò)誤裝配，因此針對(duì)long read的比對(duì)算法必須允許比對(duì)中包含空位(gap)和存在部分比對(duì)read序列［13］.為了比對(duì) long read序列，2017年，Lin等人提出的Kart［14］算法使用BWT對(duì)參考序列建立索引并進(jìn)行種子篩選，實(shí)現(xiàn)有效比對(duì)較長(zhǎng)的read序列，但當(dāng)read序列中錯(cuò)誤率與突變率較高時(shí)，Kart算法的召回率(recall rate)明顯下降.

本文采用增強(qiáng)型稀疏后綴數(shù)組(enhanced sparse suffix array)［15］對(duì)參考序列建立索引，自適應(yīng)調(diào)整種子的最小長(zhǎng)度，通過(guò)搜索參考序列與read序列之間的最大精確匹配(Maximal Exact Matches，MEMs)和超大精確匹配(Super Maximal Exact Matches，S-MEMs)作為種子查找候選區(qū)域，從而有效增加種子集合中的種子數(shù)量，以進(jìn)一步提升比對(duì)算法的召回率，實(shí)現(xiàn)映射更多的read序列.

2 改進(jìn)的long-read比對(duì)算法

2.1 Kart算法分析

Kart算法將給定read序列及其候選區(qū)域劃分為簡(jiǎn)單區(qū)域(simple pairs)和普通區(qū)域(normal pairs)，其中所有簡(jiǎn)單區(qū)域均由MEMs構(gòu)成.設(shè)read序列集R與參考序列G之間的MEMs長(zhǎng)度為L(zhǎng)k，MEMs在參考序列G中出現(xiàn)次數(shù)為Occ.Kart算法通過(guò)BWT數(shù)組搜索所有長(zhǎng)度Lk大于最小長(zhǎng)度h并且出現(xiàn)次數(shù)Occ小于50次的MEMs作為種子.當(dāng)read序列長(zhǎng)度增長(zhǎng)時(shí)，Kart算法對(duì)普通區(qū)域進(jìn)行二次篩選種子，縮短普通區(qū)域的長(zhǎng)度，以達(dá)到快速比對(duì)的目的.Kart算法中種子的最小長(zhǎng)度h的值設(shè)定在10至16之間，當(dāng)read序列長(zhǎng)度不斷增加時(shí)，設(shè)定的種子的最小長(zhǎng)度值遠(yuǎn)小于read序列的長(zhǎng)度，且通過(guò)BWT算法搜索到的MEMs之間會(huì)存在大量的重疊(overlap).另外，long read(長(zhǎng)序列)中包含許多簡(jiǎn)單和復(fù)雜的錯(cuò)誤(例如，插入刪除和結(jié)構(gòu)變異).當(dāng)long read序列中錯(cuò)誤率和突變率逐漸增加后，Kart算法的召回率明顯下降，這就意味著Kart算法能夠識(shí)別的read數(shù)量大大減少.

2.2 改進(jìn)的算法

在給出改進(jìn)算法之前，首先介紹最大精確匹配MEMs和超大精確匹配S-MEMs兩個(gè)概念［16］.

定義1.假定i，j為字符串下標(biāo)，lm為位置i至j的子串長(zhǎng)度.當(dāng)滿足如下條件時(shí)，則(i，j，lm)被稱為長(zhǎng)度n的參考序列G和長(zhǎng)度為m的read序列集R中的序列r之間的最大精確匹配(MEM):

條件2)和條件3)分別被稱為左極大性(left-maximality)和右極大性(right-maximality).如果只有條件1)和條件2)成立，則稱(i，j，lm)為左最大精確匹配.如果只有條件1)和條件3)成立，則稱(i，j，lm)為右最大精確匹配.

定義2.在reads序列與參考序列之間，不被其他最大精確匹配所包含的最大精確匹配稱為超大精確匹配(S-MEM).

為了解決Kart算法召回率較低的問(wèn)題，本文提出一種改進(jìn)的比對(duì)算法，其主要思想是:采用增強(qiáng)型稀疏后綴數(shù)組對(duì)參考序列建立索引，搜索參考序列與read序列之間的MEMs和S-MEMs，將它們作為種子定位候選區(qū)域，然后對(duì)候選區(qū)域之間的普通區(qū)域進(jìn)行無(wú)空位/空位(ungap/gap)比對(duì)，以產(chǎn)生比對(duì)結(jié)果.因?yàn)榉N子長(zhǎng)度相對(duì)于read長(zhǎng)度來(lái)說(shuō)是非常短的，所以種子的最小長(zhǎng)度h是影響比對(duì)召回率的關(guān)鍵因素之一.為了避免過(guò)多MEMs產(chǎn)生重疊，算法初始時(shí)依據(jù)參考序列長(zhǎng)度設(shè)定種子的最小長(zhǎng)度h.當(dāng)搜索MEMs和S-MEMs時(shí)，依據(jù)read的長(zhǎng)度自適應(yīng)調(diào)整種子的最小長(zhǎng)度h，獲取長(zhǎng)度大于h的MEMs和S-MEMs.從S-MEMs的定義可知，S-MEMs不與其他MEMs重疊.因此，將所有S-MEMs優(yōu)先選為種子可以減少M(fèi)EMs之間的大量重疊.特別地，當(dāng)read長(zhǎng)度越長(zhǎng)時(shí)，能夠找到的S-MEMs越多，種子之間的互相重疊區(qū)域越少.相對(duì)于MEMs來(lái)說(shuō)，S-MEMs的數(shù)量較少，若只使用S-MEMs作為種子，算法的召回率將會(huì)降低.因此，本文結(jié)合使用MEMs與S-MEMs作為種子進(jìn)行擴(kuò)展比對(duì)以提高算法召回率.

算法1描述了本文改進(jìn)的long-read比對(duì)算法(稱為sufKart算法)，算法2給出了搜索MEMs和S-MEMs作為種子的算法.

設(shè)read序列集R共有l(wèi)n條序列，每條long read序列的長(zhǎng)度為m.算法2是sufKart算法中最耗時(shí)的部分，它對(duì)每條long read搜索MEMs，然后在MEMs集合中搜索S-MEMs，搜索MEMs和S-MEMs所需的時(shí)間為O(m2+o)，其中o為右最大精確匹配的數(shù)量.算法1https://github.com/lh3/wgsim中步2循環(huán)迭代執(zhí)行l(wèi)n次，步7所需時(shí)間為O(m2+o)，步9和步10所需時(shí)間分別為t1和t2，t1和t2均小于O(m2+o).算法1所需時(shí)間為lnO(m2+o)=O(ln×m2+ln×o).改進(jìn)的 long-read比對(duì)算法 sufKart的時(shí)間復(fù)雜度為O(ln×m2+ln×o).

sufKart算法采用增強(qiáng)型稀疏后綴數(shù)組對(duì)參考序列建立索引，通過(guò)搜索參考序列與read序列之間的MEMs和SMEMs的方法，將MEMs和S-MEMs作為種子定位候選區(qū)域.在搜索MEMs和S-MEMs之前，根據(jù)read序列的長(zhǎng)度動(dòng)態(tài)調(diào)整種子的最小長(zhǎng)度，使得suf Kart算法能夠比對(duì)不同長(zhǎng)度的read序列.同時(shí)，結(jié)合使用MEMs和S-MEMs作為種子進(jìn)行擴(kuò)展的方法能夠有效減少種子之間大量的重疊，從而更準(zhǔn)確定位到相似區(qū)域.因此，sufKart算法可以明顯提升召回率.

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

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

實(shí)驗(yàn)使用的計(jì)算機(jī)配置2個(gè)10核Intel Xeon(R)CPU E5-2660，2.2GHz處理器，內(nèi)存容量128GB，操作系統(tǒng)為L(zhǎng)inux ubuntu 16.04 LTS，采用C++語(yǔ)言編程實(shí)現(xiàn)算法.

本文采用模擬數(shù)據(jù)與真實(shí)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)測(cè)試，參考序列為人類基因組(hg19).模擬數(shù)據(jù)采用Wgsim1https://github.com/lh3/wgsim模擬器自動(dòng)生成，長(zhǎng)度分別為 1000bp、1500bp、2000bp、2500bp、3000bp、4000bp和5000bp的read序列;錯(cuò)誤率分別為2%，5%，10%和15%，突變率分別為2%，5%，8%和10%.真實(shí)數(shù)據(jù)采用Ion Torrent測(cè)序平臺(tái)產(chǎn)生的reads平均長(zhǎng)度為172bp的數(shù)據(jù)集SRR946843，454測(cè)序平臺(tái)產(chǎn)生的reads平均長(zhǎng)度為574bp的數(shù)據(jù)集SRR003161和PacBio平臺(tái)產(chǎn)生的reads平均長(zhǎng)度為7116bp的數(shù)據(jù)集M130929.

在模擬數(shù)據(jù)實(shí)驗(yàn)中，本文與文獻(xiàn)［14，16］一樣使用精確度(precision rate)和召回率作為評(píng)估比對(duì)算法的性能指標(biāo).設(shè)read數(shù)據(jù)集中有N條序列，其中有M條序列被映射到參考序列(即這M條read序列在參考序列中出現(xiàn))，當(dāng)一條read的映射位置與其模擬映射位置的差值在30bp以內(nèi)，則認(rèn)定這條read被正確映射(true mapped)［14］.記 reads數(shù)據(jù)集中被正確映射到參考序列的read數(shù)量為TM，序列比對(duì)算法的精確度((precision rate)和召回率(recall rate)的計(jì)算方法如下［14］:

3.2 模擬數(shù)據(jù)實(shí)驗(yàn)結(jié)果

首先，對(duì)長(zhǎng)度分別為 1000bp、1500bp、2000bp、2500bp、3000bp、4000bp和5000bp的模擬數(shù)據(jù)read序列，對(duì)于錯(cuò)誤率與突變率不同取值的情況，實(shí)驗(yàn)測(cè)試sufKart和Kart算法的精確度與召回率.對(duì)于不同長(zhǎng)度的reads，當(dāng)錯(cuò)誤率分別等于2%、5%、10%和15%(突變率取Wgsim的默認(rèn)值)時(shí)，圖1和圖2分別給出了Kart與sufKart算法的精確度和召回率.

根據(jù)以上公式可知，如果每條read都被正確映射，那么精確度等于召回率.

在真實(shí)數(shù)據(jù)實(shí)驗(yàn)中，因?yàn)闊o(wú)法預(yù)知read的實(shí)際映射位置，且最優(yōu)的比對(duì)算法應(yīng)可以識(shí)別最大數(shù)量的堿基，所以采用運(yùn)行時(shí)間、敏感度和平均比對(duì)長(zhǎng)度作為評(píng)估指標(biāo).比對(duì)算法的敏感度(sensitivity)的計(jì)算方法如下［14］:

圖1 不同錯(cuò)誤率的模擬數(shù)據(jù)上算法的精確度Fig.1 Accuracy of algorithms for simulated reads with different error rates

從圖1和圖2中可以看出:當(dāng)錯(cuò)誤率為5%、read長(zhǎng)度分別等于1500 bp、2000bp和3000bp時(shí)，sufKart算法的精確度高于Kart算法.當(dāng)錯(cuò)誤率為5%、10%和15%時(shí)，sufKart算法的召回率整體上均高于Kart算法;特別地，當(dāng)錯(cuò)誤率分別為15%、read長(zhǎng)度為 1000bp時(shí)，sufKart算法的召回率為96.47%，而Kart算法的召回率僅為70.68%.這表明，當(dāng)read的錯(cuò)誤率逐漸增大時(shí)，本文的sufKart算法的召回率比Kart算法更高，且能映射更多的reads序列.

對(duì)于不同長(zhǎng)度的reads，當(dāng)突變率分別等于2%、5%、8%和10%(錯(cuò)誤率取Wgsim的默認(rèn)值)時(shí)，圖3和圖4分別給出了Kart與sufKart算法的精確度和召回率.

圖2 不同錯(cuò)誤率的模擬數(shù)據(jù)上算法的召回率Fig.2 Recall rate of algorithms for simulated reads with different error rates

由圖3和圖4可知:當(dāng)突變率高達(dá)8%以上時(shí)，不論是精確度還是召回率，sufKart算法幾乎均高于Kart算法.

圖3 不同突變率的模擬數(shù)據(jù)上算法的精確度Fig.3 Accuracy of algorithms for simulated reads with different mutation rates

圖4 不同突變率的模擬數(shù)據(jù)上算法的召回率Fig.4 Recall rate of algorithms for simulated reads with different mutation rates

在實(shí)際應(yīng)用中，各類測(cè)序平臺(tái)產(chǎn)生的reads不僅僅只包含簡(jiǎn)單的插入刪除錯(cuò)誤而且還包括比較復(fù)雜的結(jié)構(gòu)變異.因此，本文生成不同錯(cuò)誤率和突變率的模擬reads數(shù)據(jù)集進(jìn)行實(shí)驗(yàn).表1給出了長(zhǎng)度1000bp，錯(cuò)誤率分別為2%、5%、10%和15%，突變率為2%、5%、8%和10%時(shí)，Kart和 sufKart算法對(duì)模擬數(shù)據(jù)reads進(jìn)行實(shí)驗(yàn)獲得的精確度和召回率，其中“E02-R02”表示錯(cuò)誤率為2%，突變率為2%的reads數(shù)據(jù)集，其他類推.

從表1可以看出:當(dāng)每組模擬數(shù)據(jù)reads的(錯(cuò)誤率，突變率)分別為(2%，5%)、(2%，8%)、(5%，2%)和(5%，5%)時(shí)，sufKart算法的精確度高于Kart算法，且召回率幾乎也高于Kart算法.當(dāng)錯(cuò)誤率為10%和15%時(shí)，雖然Kart算法的精確度稍高于sufKart算法的精確度，但是sufKart算法的召回率遠(yuǎn)高于Kart算法;特別地，當(dāng)錯(cuò)誤率高達(dá)15%時(shí)，sufKart算法的召回率比Kart算法高2至3倍.此外，sufKart算法的精確度的值與召回率的值彼此更為接近，也就是說(shuō)sufKart算法能夠正確映射更多的reads序列.

表1 不同錯(cuò)誤率和突變率的模擬數(shù)據(jù)上算法的精確度和召回率Table 1 Accuracy and recall rate of algorithms with different error rates and mutation rates on simulation data

3.3 真實(shí)數(shù)據(jù)實(shí)驗(yàn)結(jié)果

本文采用3組真實(shí)數(shù)據(jù)集 SRR946843、SRR003161和M130929對(duì) sufKart和 Kart算法進(jìn)行實(shí)驗(yàn)測(cè)試.其中，SRR946843是Ion Torrent測(cè)序平臺(tái)產(chǎn)生的平均長(zhǎng)度為172bp的reads數(shù)據(jù)集，SRR003161是454測(cè)序平臺(tái)產(chǎn)生的平均長(zhǎng)度為574bp的reads數(shù)據(jù)集，M130929是PacBio測(cè)序平臺(tái)產(chǎn)生的平均長(zhǎng)度為7116bp的reads數(shù)據(jù)集.兩種算法的平均比對(duì)長(zhǎng)度和敏感度的結(jié)果如表2所示.

表2的實(shí)驗(yàn)結(jié)果表明，對(duì)于真實(shí)數(shù)據(jù)集SRR946843和SRR003161，sufKart算法在平均比對(duì)長(zhǎng)度和敏感度兩方面明顯高于Kart算法.對(duì)于M130929數(shù)據(jù)集，Kart算法的敏感度與suf Kart算法的敏感度幾乎相同;當(dāng)read序列長(zhǎng)度達(dá)到7000bp以上時(shí)，Kart算法的平均比對(duì)長(zhǎng)度大于sufKart算法的平均比對(duì)長(zhǎng)度，這是因?yàn)榇藭r(shí)sufKart算法中種子集合內(nèi)的MEMs數(shù)量不夠多，導(dǎo)致候選區(qū)域數(shù)量較少，從而降低了sufKart算法的平均比對(duì)長(zhǎng)度.在運(yùn)行時(shí)間方面，sufKart算法所需的運(yùn)行時(shí)間稍高于Kart算法.這是因?yàn)閟ufKart算法需要花費(fèi)一些時(shí)間，用于搜索reads和參考序列之間的最大精確匹配(MEMs)和超大精確匹配(S-MEMs)作為種子查找候選區(qū)域，以有效增加種子集合中的種子數(shù)量，從而使得sufKart算法比對(duì)結(jié)果的敏感度較之Kart算法有所提升.

表2 真實(shí)數(shù)據(jù)上算法的平均比對(duì)長(zhǎng)度、敏感度和運(yùn)行時(shí)間Table 2 Average alignment length，sensitivity and runtinme for algorithms on real data

綜合模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)的實(shí)驗(yàn)結(jié)果表明，本文算法sufKart的比對(duì)結(jié)果在召回率和敏感度總體上優(yōu)于Kart算法，且sufKart能夠識(shí)別更多的reads.

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

本文給出的sufKart算法的新穎之處是:采用增強(qiáng)型稀疏后綴數(shù)組對(duì)參考序列建立索引，在篩選種子階段，種子的最小長(zhǎng)度依據(jù)read長(zhǎng)度和錯(cuò)誤率進(jìn)行自動(dòng)調(diào)整，同時(shí)搜索MEMs和S-MEMs，產(chǎn)生更多有效的種子，這些種子之間重疊較少，減少了大量去除種子之間重疊的工作，能夠更準(zhǔn)確的定位候選區(qū)域.sufKart算法在篩選種子階段需要花費(fèi)較多時(shí)間搜索MEMs和S-MEMs以產(chǎn)生更多的有效種子.對(duì)于計(jì)算密集型的任務(wù)，可以采用GPU計(jì)算來(lái)加速.我們下一步的工作將借鑒文獻(xiàn)［17］的并行比對(duì)處理思想，研究設(shè)計(jì)GPU并行算法來(lái)提高種子篩選的效率，以減少long-read比對(duì)算法所需的時(shí)間.