基于高通量測序數據的基因組從頭組裝方法
【專利摘要】本發明提供了基于高通量測序數據的基因組從頭組裝方法，包括步驟：1)依據高通量測序數據構建de Bruijn圖，基于糾錯后的de Bruijn圖進行測序數據糾錯和super read組裝；2)利用super read進行初級contigs組裝；3)調取特定局部的初級contigs和reads，局部組裝，將所有的局部組裝結果合并；4)通過子圖分割算法和模擬退火算法對contigs進行排序得到最終的scaffolds。本發明通過de Bruijn圖糾錯消除高通量測序帶來的錯誤，提高了數據準確性；采用構建super read的方法提高測序讀長，顯著提升contigs長度；通過局部組裝大大提升了重復序列的處理能力。
【專利說明】基于高通量測序數據的基因組從頭組裝方法
【技術領域】
[0001] 本發明涉及一種基因組組裝方法，特別是涉及一種基于短序列測序片段的基因組 從頭組裝方法。
【背景技術】
[0002] 隨著第二代測序技術的迅速發展，測序費用的急速下降，從頭基因組測序愈發受 到研究者的青睞。但是，利用大量的短的read數據重新恢復出基因組原貌也面臨著巨大 的挑戰，而其中最為關鍵的一步就是contigs組裝。DeBruijn圖構建是圖論組裝算法的 核心，它是現在主流從頭組裝算法的核心，它是基于kmer的重疊信息來構建歐拉圖，它是 contigs構建的基石，因此本發明的開發也將基于DeBruijn圖。
[0003] 目前的contigs組裝算法都只進行一次DeBruijn圖構建，同時針對圖中的kmer 大小也是相對固定的，雖然存在一些多kmer組裝算法，但它們也都只進行一次構圖，再進 行合并。對于組裝中使用的短序列一般的組裝軟件也只是進行簡單的過濾也糾錯處理，并 不會對這些最原始的短序列進行二次加工，這也就在很大程度上限制的DeBruijn圖構建 中kmer大小的上限。因此對于不進行短序列加工的基因組組裝方法，kmer大小都比較小， 在DeBruijn圖構建中會產生較多的分支，極大的提高DeBruijn圖的復雜度，從而降低組 裝效果。
[0004] 另外，動植物基因組的一大特點就是重復序列比例較高，而重復序列會讓基因組 組裝過程中產生大量的可選位點和分支，進而提高組裝難度。目前主要有兩種主流的策略 來處理其中的部分情況：一種策略是利用大片段文庫跨過重復序列，并估計重復序列區域 大小，然后選取一個合適長度的重復序列路徑；另一種則是先回避重復序列區域，在完成初 步組裝后再回過頭來進行重序區域的組裝。從策略上講，第二種方法對于復雜基因組來說 更有效，因為它把全局問題進行了局部化，大大降低了組裝的難度。
【發明內容】
[0005] 針對現有技術存在的不足，本發明的目的是提供一種基于高通量測序數據的基因 組從頭組裝方法--GN0V0,該技術首先通過數據糾錯來處理高通量測序固有的測序錯誤， 同時通superread組裝將較短的read組裝為具有更大讀長的superread,從而部分克服 測序讀長過短的問題。其次，通過局部組裝，將全基因組上的重復序列轉變為局部的單拷貝 序列，從而大大降低了重復序列處理的難度，提高了contigs組裝的長度。
[0006] 為了實現本發明目的，本發明的一種基于高通量測序數據的基因組從頭組裝方 法--GN0V0,主要步驟為：
[0007] 1)通過高通量測序數據構建deBruijn圖（使用較小的kmer)，并進行圖糾錯處 理，并基于糾錯后的deBruijn圖進行測序數據糾錯，糾錯原理見圖4 ;
[0008] 2)基于糾錯后的deBruijn圖進行superread組裝；
[0009]3)用superread重新構建deBruijn圖（使用較大的kmer)，并進行圖糾錯處理， 對糾錯后的deBruijn圖進行拆分，得到初級contigs ;
[0010] 4)根據mate-pair的連接信息調取特定局部的初級contigs,并依據測序數據的 比對信息收集局部的reads進行局部組裝，將所有的局部組裝結果合并到一起，并進行糾 錯后拆分處理，從而得到contigs;
[0011] 5)根據mate-pair的連接信息構建scaffold連接圖，通過子圖分割算法 對contigs進行分割，并采用模擬退火算法在局部對contigs進行排序得到最終的 scaffolds。
[0012] GN0V0組裝原理流程圖見圖1。
[0013] 步驟1-4中都是以deBruijn圖為核心結構的，在GN0V0中，deBruijn圖是以哈 希的數據結構形式存在的，其構建算法為：
[0014] 1)根據基因組大小和kmer大小對哈希表進行空間分配與初始化；
[0015] 2)迭代讀取每條read，并進行編號，編號從0開始。
[0016] 3)從5'到3'端依次提取所有的kmer，并將其存C到哈希表中。如果kmer已經 存在，則只需存貯kmer的路徑信息就可以了，即存貯其前驅與后驅。如果kmer不存在，則 需要新建kmer節點，同時還需存路徑信息。
[0017] 4)存read中第一個kmer信息時，如果其在哈希表中不存在，則說明其真前驅 kmer節點不存在，至少到當前為止，它是不存在的。因此，這個時候就需要新建一個未端測 序突起節點，用于取代真實前驅kmer節點，作為該kmer節點的回溯前驅節點。
[0018] 5)在存C非第一個kmer節點時，如果發現該kmer已經存在了，并且該kmer節點 的回溯前驅節點為未端測序突起節點，則需要將該未端測序突起節點去掉，同時將該kmer 節點回溯前驅節點設置為前一個kmer節點。因為，在當前read中，該kmer不是第一個 kmer,所以其一定有一個真實前驅，即它的前一個kmer,因此，可以用真實前驅kmer節點來 代替未端測序突起節點，從而減少未端測序突起數量，進而節省部分內存。
[0019] deBruijn圖作為核心數據結構，其準確性是十分重要的，因此，GN0V0中開發了一 系列的圖糾錯處理，主要步驟包裝：l)deBruijn圖簡化處理；2)未端測序突起刪除處理； 3)泡狀路徑合并操作；4)低覆蓋度邊清除處理。
[0020] l)de Bruijn圖簡化處理：依據哈希表,對每個kmer節點進行遍歷。對于當前kmer 節點，依據其真實前驅與后驅進行延伸，如果當前kmer節點的互補節點也存在，則需要同 時依據這兩個節點進行延伸。延伸方法：沿著出邊與入邊的方向進行延伸，即沿真實前驅 與后驅進行延伸。單個方向上的延伸條件：延伸處的kmer節點（包括其互補kmer,如果存 在）有且僅有一個真實前驅，同時有且僅有一個真實后驅。單個方向上的延伸終止：對于后 驅延伸來說，延伸到的當前kmer節點有兩個或多個真實前驅，或者說有兩個或多個真實后 驅，或者說沒有真實后驅了，或者說延伸到的當前kmer節點已經在De Bruijn圖中存在了。 對于前驅延伸來說，延伸到的當前kmer節點有兩個或多個真實前驅，或者說有兩個或多個 真實后驅，或者說沒有真實前驅了，或者說延伸到的當前kmer節點已經在De Bruijn圖中 存在了。
[0021] 2)未端測序突起主要是由于read末端的測序錯誤產生的，GN0V0中未端測序突起 錯誤的判斷標準為：a)長度小于2K(K為kmer的長度）；b)必須存在高覆蓋度的等位入邊 或出邊。
[0022] 3)泡狀路徑是指由具有相同起點和終點的兩條不同的路徑構成的圖形結構，除起 點和終點外，圖形內部不存在其它任何的交叉節點。泡狀路徑主要是由雜合位點與read中 部的測序錯誤產生的，GN0V0中泡狀路徑的定義為：1)路徑長度均小于200bp ;2)路徑的相 似度大于0.8 ;3)至少有一條路徑的覆蓋度低于某個特定的閾值。泡狀路徑搜索算法的核 心算法是"Dijkstra-like breadth-first search"（Dijkstra算法是最短路徑搜索算法中 最著名的算法，"breadth-first search"表示廣度優先遍歷）。
[0023] 4)低覆蓋度邊主要是由read測序錯誤產生的，其主要的判別標準：1)覆蓋度小于 某個特定的閾值；2)邊兩端的節點均存在除當前邊外的至少一個真實前驅和至少一個真 實后驅。覆蓋度閾值的選取，對于單倍體來說，一般默認選取邊覆蓋度的均值或者說中位數 的1/2,對于雙倍體基因組來說，默認選取邊覆蓋度的均值或者說中位數的1/4。但是最好 的方法是根據覆蓋度的總體分布進行閾值的選取。
[0024] Super Read是指一條較長的序列，它是通過補齊paired-end之間的缺口或者說 是通過重疊信息連接paired-end兩端而得到的一條序列，Super Read的構建原理見圖5。 由于它是基于paired-end獲得的，因此它的長度的期望值將為文庫片段大小。由于super read是連接了雙端的read與中間的缺口，因此其長度一般比read的讀長長很多，以super read作為組裝起點具有非常大的優勢。Super read的組裝是采用深度優先算法進行路徑 搜索得到的。
[0025] 在很多分析中，都是基于單拷貝節點出發的，原因主要有：1)先從單拷貝節點出 發，組裝分析比較容易，而且出錯的概率會較小。2)有了單拷貝節點的信息，則在后期處理 重復序列時可以借助它作為基點，解決一部分重復序列組裝。
[0026] 這里假設有一條邊，它的長度為n, Xi表示以邊上的位點i為read起始位點的 read數目（注意這里邊長的實際長度為n-k+1，因為邊是以kmer為基礎的，因此i的最大 值為n-k+1)。這里假設Xi為獨立的隨機變量，它是服從期望為P的泊松分布，它的期望 P由邊的覆蓋度的分布來確定（這里是指所有邊的覆蓋度的分布情況，即總體分布）。
[0027] 根據中心極限理論，一個長度為n的邊上的Xi的期望值應該服從均值為P，標準 差為V>Z(n-l<+l)的正態分布。如果某條邊是單拷貝邊，那么Xi的平均值與P的差異就 不應太大。這里取下面的比值作為邊唯一性的判定準確：
[0028]
【權利要求】
1. 一種基于高通量測序數據的基因組從頭組裝方法，包括以下步驟： (1) 通過高通量測序數據構建de Bruijn圖，并基于糾錯后的deBruijn圖進行高通量 測序數據糾錯； (2) 采用路徑搜索算法構建super read ; (3) 利用super read重新構建de Bruijn圖，通過de Bruijn圖糾錯和拆分處理，得到 初級 conti gs ; (4) 根據pair-end和mate-pair的比對信息調取特定局部的初級contigs和reads, 進行局部組裝； (5) 利用mate-pair連接信息構建scaffold連接圖，通過子圖分割算法和模擬退火算 法對contigs進行排序得到最終的scaffolds。
2. 根據權利要求1所述的方法，其特征在于，步驟（1)中構建de Bruijn圖的kmer長 度根據覆蓋度信息進行浮動，長度在17-37之間。
3. 根據權利要求1所述的方法，其特征在于，步驟（1)所述的de Bruijn圖是以哈希的 數據結構形式存在的，其構建算法為： 1) 根據基因組大小和kmer大小對哈希表進行空間分配與初始化； 2) 迭代讀取每條read，并進行編號，編號從0開始； 3) 從5'到3'端依次提取所有的kmer，并將其存貯到哈希表中，如果kmer已經存在， 則只需存C kmer的路徑信息，即存其前驅與后驅；如果kmer不存在，則需要新建kmer節 點，同時還需存It路徑?目息； 4) 存C read中第一個kmer信息時，如果其在哈希表中不存在，則說明其真前驅kmer 節點不存在，需要新建一個未端測序突起節點，用于取代真實前驅kmer節點，作為該kmer 節點的回溯前驅節點； 5) 在存C非第一個kmer節點時，如果發現該kmer已經存在，并且該kmer節點的回溯 前驅節點為未端測序突起節點，則需要將該未端測序突起節點去掉，同時將該kmer節點回 溯前驅節點設置為前一個kmer節點。
4. 根據權利要求1所述的方法，其特征在于，步驟（1)和步驟（3)所述的de Bruijn圖 是通過以下步驟來糾錯的：l)de Bruijn圖簡化處理；2)未端測序突起刪除處理；3)泡狀路 徑合并操作；4)低覆蓋度邊清除處理。
5. 根據權利要求1所述的方法，其特征在于，步驟（2)的路徑搜索算法為深度優先路徑 搜索算法。
6. 根據權利要求1所述的方法，其特征在于，步驟（3)是使用較大的kmer構建de Bruijn 圖，kmer 大小為 75-155。
7. 根據權利要求1所述的方法，其特征在于，步驟（4)的局部組裝步驟為： 1) 將初級contigs和reads做比對，通過reads的比對結果，得到初級contigs之間的 距離信息，以及reads和初級contigs的關系，將初級contigs和reads信息讀入內存； 2) 過濾掉多拷貝（拷貝數>2)或長度較短的初級contigs,，根據初級contigs之間 的距離關系，構建scaffold連接圖，并在其中選擇相距較遠且較長的初級contigs作為種 子，得到種子后選出在種子附近一定范圍內的初級contigs ; 3) 對每個局部的初級contigs,根據比對結果選擇只有一端在初級contigs上的測序 片段，同時將處在缺口處并且測序片段覆蓋度大于0. 9的super read也選取出來； 4) 在局部構建de Bruijn圖進行局部組裝； 5) 把每個局部圖內的局部組裝結果進行合并，得到全局的組裝結果，然后進行簡化與 圖糾錯處理，從而得到最終的contigs。
8. 根據權利要求1-6任一項所述的方法，其特征在于，步驟5通過子圖分割算法和模擬 退火算法對contigs進行排序時，根據子圖的大小不同選取直接排序與模擬退火排序，子 圖< 8采用直接排序，> 8則采用模擬退火排序。
9. 根據權利要求8所述的方法，其特征在于，所述子圖分割算法包括以下步驟： 1) 依次遍歷每個contig ; 2) 檢測每個contig的長度和邊連接； 3) 對定義為邊界contig進行刪除，即進行圖拆分； 所述直接排序算法包括以下步驟： 1) 對所有的排序可能進行窮舉； 2) 選取沖突邊最少的排序作為最優排序； 模擬退火算法對contigs進行排序包括以下步驟： 1) 隨機生成一個contigs排序； 2) 根據當前溫度和隨機改變生成的新排序計算一個接收概率； 3) 生成一個隨機概率，如果該隨機概率小于接收概率則將新的contigs排序作為當前 排序； 4) 以特定步長進行降溫，同時迭代步驟2和3 ; 5) 從所有的當前排序中選取沖突邊最少的排序作為最優排序。
10. 根據權利要求1-4任一項所述的方法，其特征在于，高通量測序數據糾錯、局部組 裝過程中的路徑搜索方式為深度優先算法，同時對于圖糾錯中的路徑搜索方式需要對邊進 行加權。
【文檔編號】G06F19/18GK104239750SQ201410421844
【公開日】2014年12月24日 申請日期:2014年8月25日 優先權日:2014年8月25日
【發明者】鄭洪坤, 劉敏 申請人:北京百邁客生物科技有限公司