基于de Bruijn圖的基因組索引結(jié)構(gòu)設(shè)計(jì)

國(guó)宏哲， 王亞?wèn)|

(哈爾濱工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院， 哈爾濱 150001)

0 引 言

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

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

1 基于de Bruijn圖的索引結(jié)構(gòu)分析

1.1 基于哈希表結(jié)構(gòu)的數(shù)據(jù)存儲(chǔ)方法

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

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

Fig.1Anexampleofgenomick-mersstorageinthehashtablestructure

1.2 基于哈希結(jié)構(gòu)的序列比對(duì)方法

k-mer長(zhǎng)度越長(zhǎng)，在比對(duì)過(guò)程中需要搜索k-mer的次數(shù)越少，且速度也越快。但在實(shí)際應(yīng)用中需要嚴(yán)格控制k-mer的長(zhǎng)度，因?yàn)榈刂房臻g是隨著其序列長(zhǎng)度的增加而呈指數(shù)增長(zhǎng)。而且，k-mer長(zhǎng)度過(guò)大會(huì)使比對(duì)時(shí)seed過(guò)長(zhǎng)而較大概率地跨過(guò)mismatch或indel存在的bp位置，將真實(shí)的seed漏掉而無(wú)法得到真實(shí)的位置信息。如果k-mer長(zhǎng)度過(guò)小會(huì)使整個(gè)de Bruijn圖過(guò)密、節(jié)點(diǎn)過(guò)多，且每個(gè)節(jié)點(diǎn)對(duì)應(yīng)的入邊和出邊也較多，這會(huì)使圖上的搜索過(guò)慢。故k-mer長(zhǎng)度的選取成為一個(gè)關(guān)鍵問(wèn)題。針對(duì)此問(wèn)題，需要計(jì)算不同k-mer長(zhǎng)度在不同物種基因組上的分布。

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

結(jié)合Hash實(shí)現(xiàn)過(guò)程中的實(shí)際情況，針對(duì)人類(lèi)基因組數(shù)據(jù)使用，本文將序列映射系統(tǒng)索引和比對(duì)過(guò)程中k-mer長(zhǎng)度范圍定為21～28，將前14 bp作為地址，并將剩余bp構(gòu)建連續(xù)的數(shù)組結(jié)構(gòu)(地址空間約228* 4 B = 1 GB)。

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

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

圖2 基于hash結(jié)構(gòu)存儲(chǔ)的read比對(duì)過(guò)程

Fig.2Anillustrationofthereadalignmentbasedonthehashstructure

圖3 種子對(duì)應(yīng)位置集合合并過(guò)程

Fig.3Anillustrationoftheprocessofseedsmergingbythepositioncollections

1.3 k-mer的3層存儲(chǔ)結(jié)構(gòu)分析

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

為實(shí)現(xiàn)該索引結(jié)構(gòu)，系統(tǒng)需要遍歷一遍基因組序列。且第二層序列內(nèi)容存儲(chǔ)和第三層位置信息需要有序存儲(chǔ)，故需要先將基因組上所有k-mer進(jìn)行排序。而由于基因組序列較長(zhǎng)無(wú)法做到將所有k-mer統(tǒng)一提取排序，故采用分段排序的辦法、即按照字母表邏輯大小順序(A

圖4 索引存儲(chǔ)三層結(jié)構(gòu)

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

在比對(duì)部分首先導(dǎo)入索引文件。某一k-mer的搜索在第一層地址空間中時(shí)間是O(1)；在第二層后半部分，k-mer序列搜索過(guò)程時(shí)間是O(n)；第三層返回的位置信息數(shù)組是由第二層結(jié)構(gòu)中的指針加以控制。其做到了最快速地返回某一k-mer的所有相關(guān)信息并實(shí)現(xiàn)了位置信息的融合且有序。基于此結(jié)構(gòu)的read的seeding過(guò)程可闡釋如下。

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

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

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

2 索引結(jié)構(gòu)構(gòu)建方法

2.1 索引的基本組織結(jié)構(gòu)分析

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

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

2.2 相似種子合并操作方法

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

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

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

3 結(jié)果分析

為測(cè)試基于索引結(jié)構(gòu)的seeding(種子獲取過(guò)程)效果，研究分別統(tǒng)計(jì)模擬數(shù)據(jù)和真實(shí)測(cè)序數(shù)據(jù)上的seeding中各過(guò)程產(chǎn)生的seeds數(shù)，并使用BWA比對(duì)軟件中的基于BWT索引結(jié)構(gòu)的索引和基于BWT MEM索引結(jié)構(gòu)(maximal exact matches，MEM)的索引進(jìn)行比較分析。實(shí)驗(yàn)中擬使用Mason模擬器分別模擬生成基于Illumina平臺(tái)的長(zhǎng)度為100 bp 和 250 bp的reads的數(shù)據(jù)集(Sim-i100和Sim-i250)。

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

表1模擬數(shù)據(jù)集上的索引結(jié)構(gòu)seeding統(tǒng)計(jì)結(jié)果分析

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 測(cè)序數(shù)據(jù)集上的索引結(jié)構(gòu)seeding統(tǒng)計(jì)結(jié)果

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

研究使用來(lái)源于千人基因組計(jì)劃中的樣本NA12878的Illumina測(cè)序平臺(tái)生成的ERR174324真實(shí)數(shù)據(jù)集測(cè)試各索引結(jié)構(gòu)seeding情況。由于真實(shí)數(shù)據(jù)中的測(cè)序錯(cuò)誤，也加入了更多、更復(fù)雜變異，從表2可看出，整體上真實(shí)數(shù)據(jù)比模擬數(shù)據(jù)上產(chǎn)生更多的seeds。所有的Valuedseedsratio相對(duì)模擬數(shù)據(jù)降低，同樣由于測(cè)序錯(cuò)誤等問(wèn)題使索引生成很多錯(cuò)誤seeds。在該情況下，索引仍然表現(xiàn)出較好的seeding結(jié)果。

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

表3測(cè)序數(shù)據(jù)集上索引結(jié)構(gòu)的不同hits數(shù)量閾值下seeding統(tǒng)計(jì)結(jié)果

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 結(jié)束語(yǔ)

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

本文首次提出基于hash思想的用于k-mer存儲(chǔ)的3層索引結(jié)構(gòu)，并以此為基礎(chǔ)提出基于de Bruijn圖的索引結(jié)構(gòu)來(lái)組織基因組中的重復(fù)序列。

同時(shí)，提出該索引的基本構(gòu)建方法和數(shù)據(jù)的存儲(chǔ)格式。之后提出的結(jié)構(gòu)上的PRP集合合并操作、相同的unipath上的種子合并方法和unipath上的相同序列合并方法，表明索引可以極大減少用于extension過(guò)程計(jì)算的候選種子的數(shù)量。

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