增強型稀疏后綴數(shù)組索引的高錯誤率reads比對

韋 好，鐘 誠

(廣西大學計算機與電子信息學院廣西高校并行分布式計算技術重點實驗室，南寧530004)

E-mail:394178293@qq.com

1 引言

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

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

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

2 改進的long-read比對算法

2.1 Kart算法分析

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

2.2 改進的算法

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

定義1.假定i，j為字符串下標，lm為位置i至j的子串長度.當滿足如下條件時，則(i，j，lm)被稱為長度n的參考序列G和長度為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算法召回率較低的問題，本文提出一種改進的比對算法，其主要思想是:采用增強型稀疏后綴數(shù)組對參考序列建立索引，搜索參考序列與read序列之間的MEMs和S-MEMs，將它們作為種子定位候選區(qū)域，然后對候選區(qū)域之間的普通區(qū)域進行無空位/空位(ungap/gap)比對，以產生比對結果.因為種子長度相對于read長度來說是非常短的，所以種子的最小長度h是影響比對召回率的關鍵因素之一.為了避免過多MEMs產生重疊，算法初始時依據(jù)參考序列長度設定種子的最小長度h.當搜索MEMs和S-MEMs時，依據(jù)read的長度自適應調整種子的最小長度h，獲取長度大于h的MEMs和S-MEMs.從S-MEMs的定義可知，S-MEMs不與其他MEMs重疊.因此，將所有S-MEMs優(yōu)先選為種子可以減少MEMs之間的大量重疊.特別地，當read長度越長時，能夠找到的S-MEMs越多，種子之間的互相重疊區(qū)域越少.相對于MEMs來說，S-MEMs的數(shù)量較少，若只使用S-MEMs作為種子，算法的召回率將會降低.因此，本文結合使用MEMs與S-MEMs作為種子進行擴展比對以提高算法召回率.

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

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

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

3 實驗

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

實驗使用的計算機配置2個10核Intel Xeon(R)CPU E5-2660，2.2GHz處理器，內存容量128GB，操作系統(tǒng)為Linux ubuntu 16.04 LTS，采用C++語言編程實現(xiàn)算法.

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

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

3.2 模擬數(shù)據(jù)實驗結果

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

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

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

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

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

對于不同長度的reads，當突變率分別等于2%、5%、8%和10%(錯誤率取Wgsim的默認值)時，圖3和圖4分別給出了Kart與sufKart算法的精確度和召回率.

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

由圖3和圖4可知:當突變率高達8%以上時，不論是精確度還是召回率，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

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

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

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

3.3 真實數(shù)據(jù)實驗結果

本文采用3組真實數(shù)據(jù)集 SRR946843、SRR003161和M130929對 sufKart和 Kart算法進行實驗測試.其中，SRR946843是Ion Torrent測序平臺產生的平均長度為172bp的reads數(shù)據(jù)集，SRR003161是454測序平臺產生的平均長度為574bp的reads數(shù)據(jù)集，M130929是PacBio測序平臺產生的平均長度為7116bp的reads數(shù)據(jù)集.兩種算法的平均比對長度和敏感度的結果如表2所示.

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

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

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

4 結束語

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