國宏哲， 王亞東

(哈爾濱工業(yè)大學 計算機科學與技術學院， 哈爾濱 150001)

0 引 言

HTS技術的一項重要挑戰(zhàn)是將測序序列快速且準確地比對到一個或多個參考基因組上。隨著HTS技術的普及而產生的數據量越來越大，序列比對是HTS分析中計算量最大的步驟之一。read比對過程關于速度的要求正日漸提升。而且隨著越來越多的新基因組的陸續(xù)面世，更新的HTS數據分析需要將reads比對到多個基因組、而不是單基因組上。一個物種的典型鏈的基因組可以作為金標準，再利用來源于宏基因組樣本的reads比對到這些參考基因組去識別物種類型和病原體類型[1-3]。同時，de Bruijn圖是代表和組織基因組序列的基本數據結構，并已廣泛應用于各類序列分析的學術領域中，如de novo基因組裝、高通量測序(NGS)序列比對、泛基因組分析、宏基因組學分類、轉錄本同種型鑒定和定量分析、NGS序列校正等。這些也是許多基因組學研究的基礎。

時下，多數位居前沿的先進比對軟件都是根據單基因組進行設計并優(yōu)化的。同時也有很多方法被提出索引pan-genome[4-5]，且都是通過共線性多序列比對技術組織基因組。大部分索引pan-genome的方法只能用于小規(guī)模的變異研究，例如軟件GenomeMapper[4]和BWBBLE[5]只能整合SNP和短Indel變異。同時，基于共線性結構的比對速度也較慢，其速度比目前較快的比對軟件慢幾十倍[5-6]。該方法的適用范圍僅限于高度相似的多基因組，但還不能針對來源于不同鏈或不同物種的多基因組發(fā)揮效用。這些非同源性的基因組中經常包含結構變異且共線性結構不能提供諸如大型結構變異的非線性序列變換。共線性結構本身也影響序列比對結果，不適合將序列比對到多個不同基因組。產生式壓縮后綴數組索引[6]會隨著變異數量增加而呈指數增長，使其很難處理在重復區(qū)域的變異信息[7]。增廣多基因圖結構(population reference graph， PRG)[8]可以更好支持不同種類的變異。該方法是將所有reads構建成de Bruijn圖和PRG進行比對。其本質是將基因組序列比對和從頭拼接相結合。該方法主要是為局部區(qū)域而設計適合解決基因組上復雜區(qū)域上的相關問題并提高了基因組復雜區(qū)域的比對質量。但目前尚未見到有效的方法索引全基因組的PRG結構，還不能將此方法延伸至全基因組上的序列比對。為解決以上方法的限制并優(yōu)化處理基因組上大量重復片段引起的問題，本文提出一個基于de Bruijn圖的索引結構DBG-index來有效組織重復序列，從而實現單基因組和多基因組的索引和序列比對計算。

1 基于de Bruijn圖的索引結構分析

1.1 基于哈希表結構的數據存儲方法

De Bruijn圖結構可以利用hash table數據結構進行存儲，即將各個k-mer節(jié)點提取指定bp數的2-bit位表示(A/C/G/T分別編碼成二進制00/01/10/11)作為地址空間中的地址值。該地址記錄k-mer對應入邊信息和出邊信息，及指向其對應的位置信息的地址，如圖1所示。利用de Bruijn圖構建索引并實現序列比對的優(yōu)勢如下。某一k-mer序列在基因組中可能出現在多個位置上，de Bruijn圖中某一節(jié)點表示的k-mer只出現一次且該節(jié)點對應的所有位置信息記錄在一個數據結構中。在基于后綴樹結構和基于BWT結構構建的索引中，每個后綴在基因組中對應唯一的一個位置，其位置信息是分散的，但de Bruijn圖中節(jié)點將位置信息集合在一起且允許有序?？梢赃M行相鄰節(jié)點對應位置集合的合并操作以快速篩選提取真實位置信息，并驗證當前路徑是否正確。通過對現在映射系統原理的分析，發(fā)現基于seed-and-extension思想方法[9-11]比對過程中對候選位置進行的extension局部比對是相對耗費時間的計算步驟。如何有效地縮小候選位置的范圍、減少extension的次數對系統速度實現至關重要。而基于de Bruijn圖結構的遍歷過程會大大減少候選extension的位置集合數量。

de Bruijn圖結構的一個特點是可以實現雙向圖，即某一個節(jié)點既可以記錄其出邊，也可以記錄其入邊(基因組上k-mer序列對應的上一個bp)。該特性可以使種子的延伸方法更加靈活，可以采用雙向延伸的方法，從而更快速獲得較長的種子和真實的候選位置。

Fig.1Anexampleofgenomick-mersstorageinthehashtablestructure

1.2 基于哈希結構的序列比對方法

k-mer長度越長，在比對過程中需要搜索k-mer的次數越少，且速度也越快。但在實際應用中需要嚴格控制k-mer的長度，因為地址空間是隨著其序列長度的增加而呈指數增長。而且，k-mer長度過大會使比對時seed過長而較大概率地跨過mismatch或indel存在的bp位置，將真實的seed漏掉而無法得到真實的位置信息。如果k-mer長度過小會使整個de Bruijn圖過密、節(jié)點過多，且每個節(jié)點對應的入邊和出邊也較多，這會使圖上的搜索過慢。故k-mer長度的選取成為一個關鍵問題。針對此問題，需要計算不同k-mer長度在不同物種基因組上的分布。

很多物種的基因組數據的k-mer長度變化趨勢，開始快速增長到一個值后趨于穩(wěn)定。其中，hg18人類基因組數據快速增長到20后，開始緩慢增長，一直到40后趨于穩(wěn)定。

結合Hash實現過程中的實際情況，針對人類基因組數據使用，本文將序列映射系統索引和比對過程中k-mer長度范圍定為21～28，將前14 bp作為地址，并將剩余bp構建連續(xù)的數組結構(地址空間約228* 4 B = 1 GB)。

基于Hash索引結構進行序列比對的基本原理如圖2所示。從read上從頭向后提取k-mer序列轉換成地址編碼，在Hash索引中進行搜索獲得對應的獲選位置，在對應位置提取序列和read做局部或全局比對。這也是基于seed-and-extension的經典流程。

基于de Bruijn圖的比對過程與此有所不同，在獲得某節(jié)點的位置信息后會有合并過程，通過此有序位置數據間的合并過程可篩去大量無意義的位置信息，如圖3所示。某一位置集合里的位置數值加上read上偏移，在相鄰節(jié)點對應位置集合中搜索找到合理的位置數值，從而實現兩集合間的合并，按照此方法依次向下合并。在合并過程中考慮到mismatch或indel存在的可能性，根據測序序列的錯誤率和各個類型變異的發(fā)生概率設定一個誤差范圍。

圖2 基于hash結構存儲的read比對過程

Fig.2Anillustrationofthereadalignmentbasedonthehashstructure

圖3 種子對應位置集合合并過程

Fig.3Anillustrationoftheprocessofseedsmergingbythepositioncollections

1.3 k-mer的3層存儲結構分析

各部分索引結構設計思想如圖4所示。其中包含一個3層結構：第一層是hash的地址空間保存k-mer的前14 bp序列信息；第二層記錄后面剩余的k-mer信息，空余出來的信息位可用于記錄該k-mer是否有反轉操作等信息(用于雙向seed延伸操作使用)及入邊和出邊信息；第三層是k-mer對應的有序位置信息數組。

為實現該索引結構，系統需要遍歷一遍基因組序列。且第二層序列內容存儲和第三層位置信息需要有序存儲，故需要先將基因組上所有k-mer進行排序。而由于基因組序列較長無法做到將所有k-mer統一提取排序，故采用分段排序的辦法、即按照字母表邏輯大小順序(A

圖4 索引存儲三層結構

需要注意在臨時文件中也記錄了k-mer對應的位置偏移信息及其對應的入邊和出邊信息。然后按順序重新讀取每個臨時文件，即某一時刻只有一個臨時文件的全部k-mer調取至內存中，對這些k-mer生成排序。在此過程中也可以同時構建索引的3層結構：每當臨時文件中k-mer排完序后從頭遍歷該有序結構，每當遇到一個不同k-mer序列，即將該k-mer添加到Hash地址空間中并添加入邊和出邊信息；若當前k-mer序列與前一個k-mer序列相同，則在k-mer序列對應的位置信息地址后添加新位置信息。

在比對部分首先導入索引文件。某一k-mer的搜索在第一層地址空間中時間是O(1)；在第二層后半部分，k-mer序列搜索過程時間是O(n)；第三層返回的位置信息數組是由第二層結構中的指針加以控制。其做到了最快速地返回某一k-mer的所有相關信息并實現了位置信息的融合且有序?；诖私Y構的read的seeding過程可闡釋如下。

(1)從read頭入手，提取k-mer作為seed開始在Hash索引結構上根據各節(jié)點位置集合做延伸停止后，從停止處的下一個bp開始再延伸。

(2)按此方式一直到read遍歷完成，將所有獲得的seed對應的位置結合再做進一步合并，檢驗是否有seed間可以合并的可能。

(3)將最終得到的seed按種子長度逐一排序，依次提取seed對應位置集合中位置在基因組上的序列做局部比對或全局比對；若某2個相鄰seed得到的比對分數收斂或設定某一固定的結束條件，則終止局部比對的extension環(huán)節(jié)輸出比對結果。

2 索引結構構建方法

2.1 索引的基本組織結構分析

本文構建一個(Reference de Bruijn Graph)結構去組織一個或多個reference并用到基于hash table的數據索引結構DBG-index去索引de Bruijn圖上的所有unipaths，而不是索引原始基因組。DBG-index組織和索引reference示例，一個或多個reference被組織到。一個k-mer的所有拷貝被定位到de Bruijn圖上的同一個節(jié)點，基因組上序列相同的重復片段被定位到相同的unipath上。該新基因組的索引方法提供2個基本操作去識別和合并相似種子，并且利用de Bruijn圖上unipath的性質識別相同的局部序列。此處unipath定義滿足如下條件：路徑起始節(jié)點的入度大于1、且出度是1；路徑的結束節(jié)點入度是1、且出度大于1；路徑中間節(jié)點入度是1、且出度是1。索引可以同時處理多個相似種子和相同局部序列，從而有效解決基因組的重復性帶來的問題。

比對系統是利用索引相關功能設計推出了在圖索引上的基于seed-and-extension思想的序列映射算法。算法很好處理基因組上的重復片段，從而極大減少了序列比對的時間開銷。用戶可以根據自定義k-mer長度構建reference的de Bruijn圖。de Bruijn圖上的所有節(jié)點和邊是通過reference上的所有k+1-mers計算獲得。遍歷de Bruijn圖的hash table結構，若當前k-mer為一個起始節(jié)點、即開始根據其出邊獲得其相鄰的下一個k-mer，在hash table結構中定位該k-mer，并按此方法循環(huán)直至當前k-mer為結束節(jié)點，此時可形成一個unipath序列并開始下一個unipath的生成。按此方式遍歷完hash table即可生成所有的unipath。一個unipath在基因組上可能有多個拷貝，故需要記錄每一個unipath對應的所有拷貝在基因組上的真實起始位置。索引結構包括以下3個主要結構：線性表記錄所有unipath序列且為每一個unipath指定一個標識；一個hash table索引unipath上的所有k-mers；線性表記錄每個unipath的所有拷貝的起始位置。而且通過該索引結構可實現給定k-mer獲得其所在的unipath和在unipath的偏移位置；給定一個unipath和某偏移位置獲得在基因組上原始拷貝位置集合。

2.2 相似種子合并操作方法

根據de Bruijn圖的原理，某一個k-mer在圖結構上只出現一次，任意一個指定k-mer所在的unipath和坐標是唯一的。且某個k-mer在基因組上的位置集合可通過線性表聯合計算獲得。利用以上功能接口可以通過實現以下2個主要操作(相同unipath上的相似seeds合并；unipath上相同序列合并)去合并相似種子、以及減少相同的extension計算。

某一k-mer在reference上的所有拷貝位于圖結構上的同一節(jié)點，該性質表明一個k-bp長的seed在reference上的映射等同于連接seed和圖結構上的特定節(jié)點。seed的命中是一段短片段精確比對到一個unipath的多個拷貝上。通過以上數據結構可計算某拷貝在unipath上的起始位置。

某一read映射到基因組上的某一可能位置可以被估算成某unipath拷貝在基因組上位置加上對應seed在unipath上的位置偏移。多個seeds在同一個unipath則可能有2個相似的候選位置集合且可以將相似位置集合予以合并。可以利用線性表判斷任意2個seed是否映射到同一個unipath，如果是同一unipath，可以通過發(fā)現seeds間的相似性從而進行相似seeds的合并而不需要分別計算每一個命中位置。在此基礎上，考察某一seed在其unipath上的偏移和到unipath的2個端點的距離，從而判斷該unipath序列是否容納read序列進行局部比對。如果可以，unipath對應的所有拷貝和read的所有extension計算就都將歸結為read和unipath序列本身的局部比對，而不需要對seed的每一個命中分別進行extension計算。

3 結果分析

為測試基于索引結構的seeding(種子獲取過程)效果，研究分別統計模擬數據和真實測序數據上的seeding中各過程產生的seeds數，并使用BWA比對軟件中的基于BWT索引結構的索引和基于BWT MEM索引結構(maximal exact matches，MEM)的索引進行比較分析。實驗中擬使用Mason模擬器分別模擬生成基于Illumina平臺的長度為100 bp 和 250 bp的reads的數據集(Sim-i100和Sim-i250)。

模擬數據集上的索引結構seeding統計結果和測序數據集上的索引結構seeding 統計結果可分別詳見表1和表2。由表1可見，索引結構對應的seeding過程中第一部分表示初始產生的seeds個數；第二部分表示相同unipath上seeds合并后的seeds數；第三部分表示unipath上的相同序列合并后的seeds數。BWT MEM索引對應的seeding過程中第一部分表示初始產生的seeds數；第二部分表示經過MEM計算后的seeds數。Filtering過程表示各索引結構自身定義的一些seeds的過濾規(guī)則，例如，索引中將最長seed長度小于當前read長度一半的seeds過濾掉。不同的過濾規(guī)則篩選出的seeds會有所不同。模擬數據集上Valuedseedsratio定義為產生正確比對位置的seeds占經過filtering后的seeds總數的比例；測序數據集上Valuedseedsratio定義為產生合理比對位置的seeds占經過filtering后的seeds總數的比例。

表1模擬數據集上的索引結構seeding統計結果分析

Tab.1StatisticalanalysisofDBG-indexbasedseedingonsimulationdataset

DatasetMethodSeeding #/×106Filtering #/ ×106 Valued seeds ratio/%Sim-i100DBG-index4.904.604.504.1048.50BWT-index178.90178.901.10BWT MEM-index32.6019.2012.5015.90Sim-i250DBG-index11.309.908.905.7034.50BWT-index871.70871.700.23BWT MEM-index68.4044.308.1024.80

表2 測序數據集上的索引結構seeding統計結果

由表1和表2可看出，各索引結構隨著各個計算過程seeds逐漸減少。模擬數據集和真實數據集上A索引在各個過程生成的seeds最少且最后的Valuedseedsratio最高。BWT索引在各個過程生成的seeds最多且最后的Valuedseedsratio最低。因為MEM方法有種子長度擴展和篩選功能，BWT MEM索引相對BWT索引生成更少seeds。A索引的seeding的有效性是BWT索引的幾十倍。相對于BWT MEM索引生成的seeds數目，A索引同樣表現出優(yōu)勢，且在Sim-i100上效果優(yōu)勢更加明顯。因為read長度較短時可用的seeds數也相對較少，由此即說明當seeds數較少時，索引能更快速地生成正確seeds。這對于后續(xù)能夠高效實現序列映射算法應用將十分重要。

研究使用來源于千人基因組計劃中的樣本NA12878的Illumina測序平臺生成的ERR174324真實數據集測試各索引結構seeding情況。由于真實數據中的測序錯誤，也加入了更多、更復雜變異，從表2可看出，整體上真實數據比模擬數據上產生更多的seeds。所有的Valuedseedsratio相對模擬數據降低，同樣由于測序錯誤等問題使索引生成很多錯誤seeds。在該情況下，索引仍然表現出較好的seeding結果。

在seeding過程中每個seed都有其對應的命中(hits)數，表示在基因組上對應的位置個數。由于基因組上有大量重復片段，一個seed通常對應多個hits。在映射算法中考慮到速度的要求，需要對hit設定一個閾值以實現整個系統在準確度或敏感度和速度上的平衡。為此,研究將在測序數據集上測試不同hits數對seeding情況的影響。研究得出人類基因組上的相關統計結果見表3。其中，Hits-n表示hits數上限閾值是n。隨著n的增大，seeds數增大且Valuedseedsratio降低。可以觀察到在Hits-300～-700范圍內，seeds數增長平穩(wěn)，且在Hits數大于1 000時seeds數趨于平穩(wěn)。此類統計分析有利于輔助映射算法選擇默認hits閾值。

表3測序數據集上索引結構的不同hits數量閾值下seeding統計結果

Tab.3StatisticsofDBG-indexbasedseedingwithvarioushitsvaluethresholdsonsequencingdataset

statisticsHits-100Hits-300Hits-500Hits-700Hits-900Hits-1 100Seeds #/×1071.161.271.311.341.371.37Valued seeds ratio/%17.2015.7015.2014.9014.5014.50

4 結束語

本文首先闡述基于hash table結構的一般基因組索引結構和基于此結構的比對方法。接下來分析了不同長度的k-mer在基因組上的分布情況。然后探討了seed-and-extension的基本思路，并發(fā)現根據位置集合對種子進行合并可以有效減少extension的計算量。

本文首次提出基于hash思想的用于k-mer存儲的3層索引結構，并以此為基礎提出基于de Bruijn圖的索引結構來組織基因組中的重復序列。

同時，提出該索引的基本構建方法和數據的存儲格式。之后提出的結構上的PRP集合合并操作、相同的unipath上的種子合并方法和unipath上的相同序列合并方法，表明索引可以極大減少用于extension過程計算的候選種子的數量。

本文分別在模擬數據集和真實數據集上測試索引結構和基于BWT的索引結構的seeding效果。通過比較分析seeding過程中各步驟中的seeds數量情況，發(fā)現de Bruijn圖結構能減少更多的候選seeds，同時保持更多的有效seeds比率?；赿e Bruijn圖思想構建基因組索引可以作為研究精準高效的基因組序列映射算法的基礎，同時為實現多基因組上的序列比對提供研究思路。