專利名稱：一種基于重復序列識別的全基因組測序數(shù)據(jù)的拼接方法
技術(shù)領(lǐng)域：
本發(fā)明涉及一種基于重復序列識別的全基因組測序數(shù)據(jù)的拼接方法，屬基因工程技術(shù)領(lǐng)域。
本發(fā)明提出的基于重復序列識別的全基因組測序數(shù)據(jù)的拼接方法，包括以下步驟(1)設(shè)定一個最小的DNA片段長度為15bp-20bp，計算散彈法測序數(shù)據(jù)中非重復片段出現(xiàn)的概率分布
下列公式中各參數(shù)的含意為G基因組總長，L測序平均有效讀長，N成功測序反應數(shù)，F(xiàn)識別最小片段長度，定義一個隨機變量Yik描述用散彈法對全基因組測序中上述指定長度DNA片段出現(xiàn)K次的事件 若某點開始的片段出現(xiàn)次數(shù)為k，則有k個測序片段的起點應在基因組上[I-L+F，i]區(qū)間內(nèi)，而其它N-k個測序片段的起點不在此區(qū)間內(nèi)，這一區(qū)間的長度為L-F，若所有測序片段起點在基因組上隨機分布，則根據(jù)古典概型，上述隨機變量等于1的概率為P(Yik=1)=CNk(L-F+1G)k(1-L-F+1G)N-k---(1)]]>一次測序中出現(xiàn)次數(shù)為k的片段的平均個數(shù)可表示為E(Yk)=E(Σi=1GYik)=G·GNk(L-F+1G)k(1-L-F+1G)N-k----(2)]]>使用下式作為一次測序中出現(xiàn)次數(shù)為k的片段出現(xiàn)概率的估計值；Pk＝E(Yk)/G (3)(2)計算重復片段出現(xiàn)的概率分布設(shè)片段為一個有m個拷貝的重復序列，出現(xiàn)在基因組中的m個不同位置，在散彈法測序數(shù)據(jù)集中的出現(xiàn)次數(shù)是所有位置出現(xiàn)次數(shù)的和，用Gmk表示有m個拷貝的重復序列在一次測序中的出現(xiàn)次數(shù)為k的概率，則上述關(guān)系可用數(shù)學公式表示為Gm0=P0m]]>Gm1=Cm1·P1·P0m-1]]>Gm2=Cm2·P12·P0m-2+Cm1·P2·P0m-1]]>Gm3=Cm3·P13·P0m-3+Cm2·C21·P1·P2·P0m-2+Cm1·P3·P0m-1]]>Gmj+＝1-Gm0-Gm1…-Gmj-1其中Gmj+表示出現(xiàn)次數(shù)為j和更多的概率；(3)重復序列的識別選取非重復片段出現(xiàn)概率最接近0.3％左右的次數(shù)為重復片段的判別標準，超過這一標準的片段就認為它屬于重復序列，否則就是非重復序列；(4)首先屏蔽重復序列，將上述散彈法測序數(shù)據(jù)中與識別出的片段相同的堿基改寫為N，屏蔽后剩余長度超過一個定值的測序數(shù)據(jù)仍進入拼接過程；(5)若目標基因組大小為1百萬-3千萬堿基，則屏蔽去重復序列后不分組直接進入拼接，若目標基因組明顯大于上述范圍，則需要按照測序讀數(shù)之間的關(guān)連進行分組，例如可以將參加拼接的讀數(shù)隨機分為若干組，每組讀數(shù)個數(shù)在5至10萬個之間，每組數(shù)據(jù)進行初步拼接，對拼接得到的大片段進行比較，同源性高的聚為一組，把組成它們的讀數(shù)重新拆出來放在一起，再次進行拼接；(6)將得到的大片段中的N恢復成原有堿基，并利用同一克隆正反向測序的信息找出相關(guān)的大片段以及可能出現(xiàn)在它們之間的讀數(shù)，并將其連接起來(7)所有能連接的片段都連接以后，再使用正反向測序信息把大片段排好順序，即得到目標基因組的工作框架圖。
本發(fā)明的基于重復序列識別的全基因組測序數(shù)據(jù)的拼接方法，具有提高效率、能處理復雜基因組、明顯減少錯誤拼接出現(xiàn)的概率、減少大量前期生物學實驗準備等優(yōu)點。采用本方法進行了水稻基因組的拼接工作，結(jié)果顯示本方法完全可以勝任水稻這樣復雜基因組的拼接工作，在只進行了4.2倍基因組總長測序的情況下，拼接得到的大片段已經(jīng)覆蓋了基因組中90％以上的基因，拼接錯誤率在1％左右，這大致相當于已有技術(shù)中塞萊拉公司對果蠅基因組進行13倍基因組總長測序后得到的結(jié)果，果蠅基因組的重復序列實際上明顯少于水稻，因此其拼接難度也大大小于水稻。另外，使用本發(fā)明的方法對1％人類基因組數(shù)據(jù)拼接結(jié)果表明，可以節(jié)省93％的計算機機時和84％的計算機內(nèi)存空間。
圖2表示不同測序量各種拷貝數(shù)重復序列被選出的概率比較，標準的選擇使拷貝1序列選出概率保持在0.3％左右。
圖3顯示各插入片段均取10X冗余度時對不同長度洞均可覆蓋兩次的概率。
下列公式中參數(shù)意義G基因組總長，L測序平均有效讀長N成功測序反應數(shù)，F(xiàn)識別最小片段長度。
計算霰彈法測序中非重復的小片斷的出現(xiàn)次數(shù)定義一個隨機變量Yik描述用散彈法對全基因組測序中上述指定長度DNA片段出現(xiàn)K次的事件 若某點開始的片段出現(xiàn)次數(shù)為k，則有k個測序片段的起點應在基因組上[I-L+F，i]區(qū)間內(nèi)，而其它N-k個測序片段的起點不在此區(qū)間內(nèi)，這一區(qū)間的長度為L-F，若所有測序片段起點在基因組上隨機分布，則根據(jù)古典概型，上述隨機變量等于1的概率為P(Yik=1)=CNk(L-F+1G)k(1-L-F+1G)N-k---(1)]]>一次測序中出現(xiàn)次數(shù)為k的片段的平均個數(shù)可表示為E(Yk)=E(Σi=1GYik)=G·GNk(L-F+1G)k(1-L-F+1G)N-k----(2)]]>使用下式作為一次測序中出現(xiàn)次數(shù)為k的片段出現(xiàn)概率的估計值；Pk＝E(Yk)/G(3)在實際使用重復序列識別程序時，片段長度是根據(jù)基因組大小選擇的。對于水稻基因組來說，其基因組大小大約是430Mb，小片段總數(shù)大約是108數(shù)量級，因此選擇20bp長的小片段。此時共有約1012種不同小片段，可保證不會由于隨機因素而在基因組中出現(xiàn)相同的。如果考慮的是細菌基因組，小片段總數(shù)約為106數(shù)量級，其長度可縮短為15bp，仍可保證不會由于隨機而出現(xiàn)相同。
計算重復序列中指定長度片段的出現(xiàn)次數(shù)分析推導過程不難看出上述概率都是非重復序列的概率，因為實際上假設(shè)每個片段都只出現(xiàn)在基因組的一個地方。如果片段是有m個拷貝的重復序列，它將出現(xiàn)在基因組的m個不同位置。在霰彈法測序數(shù)據(jù)集中看到的出現(xiàn)次數(shù)將是所有這些位置測序出現(xiàn)次數(shù)的和。例如測序集中出現(xiàn)次數(shù)為0，意味著所有m個位置出現(xiàn)次數(shù)都要為0；測序集中出現(xiàn)次數(shù)為1，則只有一個位置覆蓋為1，其它都要為0；而測序集中出現(xiàn)次數(shù)為2，則可能只有一個位置出現(xiàn)次數(shù)為2，其它都為0；或有兩個位置為深度1，其它均為0；等等。用Gmk表示有m個拷貝的重復序列在一次測序中的出現(xiàn)次數(shù)為k的概率，則上述關(guān)系可用數(shù)學公式表示為設(shè)片段為一個有m個拷貝的重復序列，出現(xiàn)在基因組中的m個不同位置，在散彈法測序數(shù)據(jù)集中的出現(xiàn)次數(shù)是所有位置出現(xiàn)次數(shù)的和，用Gmk表示有m個拷貝的重復序列在一次測序中的出現(xiàn)次數(shù)為k的概率，則上述關(guān)系可用數(shù)學公式表示為Gm0=P0m]]>Gm1=Cm1·P1·P0m-1]]>Gm2=Cm2·P12·P0m-2+Cm1·P2·P0m-1]]>Gm3=Cm3·P13·P0m-3+Cm2·C21·P1·P2·P0m-2+Cm1·P3·P0m-1]]>Gmj+＝1-Gm0-Gm1…-Gmj-1其中Gmj+表示出現(xiàn)次數(shù)為j和更多的概率；(4)識別重復序列根據(jù)上述概率，算出特定測序條件下(指基因組長度、總測序量、平均讀長等參數(shù))非重復序列和不同拷貝數(shù)重復序列中片段表現(xiàn)為一定測序出現(xiàn)次數(shù)的概率。然后確定一個適當?shù)某霈F(xiàn)次數(shù)標準，超過這一標準的片段就假定它屬于重復序列，否則就是非重復序列。根據(jù)各Gmk值算出不同拷貝數(shù)重復序列在這一判斷標準下漏掉的概率，如果已知基因組中各種拷貝數(shù)重復序列所占比例，則還可知道挑出的重復序列中各種拷貝數(shù)所占的比例。
為節(jié)省篇幅，只給出

圖1顯示測序覆蓋度為1X和4X時選擇不同標準各種拷貝重復序列被選出的概率。實際工作中選擇標準最重要的指標是拷貝數(shù)為1的序列(即非重復序列)被選出的概率。由于非重復序列在基因組中一般要占到總長的2/3或更多，它被選出的概率必須充分小，以便保證選出的絕大多數(shù)確實是重復序列。在不同測序覆蓋度下，本發(fā)明確定這一概率保持在0.3％左右。表1.表示不同測序量下的重復序列識別標準。圖2表示在這一標準下不同測序量各種拷貝數(shù)重復序列被選出的概率比較。從圖2中可見測序量達到4X以上時進一步增加測序量對重復序列識別的改善已不太明顯，此時拷貝數(shù)在5以上的重復序列基本上都可識別。
表1.不同測序量下的重復序列識別標準

有了上述統(tǒng)計模型，就可以根據(jù)3-4X測序數(shù)據(jù)識別出絕大部分拷貝數(shù)在5以上的重復序列。這使我們可以建立一套測序流程，處理高等生物復雜基因組全基因組鳥槍法測序數(shù)據(jù)。該流程主要步驟如圖3所示。
屏蔽重復序列的方法是把測序數(shù)據(jù)中與識別出的片段相同的堿基改寫為N。屏蔽后剩余長度超過100bp的測序數(shù)據(jù)仍進入拼接過程。屏蔽掉重復序列后大大減低了拼接的復雜性，使phrap等常用軟件可以處理更大的數(shù)據(jù)集。如果目標基因組大小在數(shù)百萬到數(shù)千萬堿基的范圍，屏蔽掉重復序列后可不分組直接投入拼接。但若目標基因組有數(shù)億或更多堿基，則需要按照測序讀數(shù)之間的關(guān)連進行分組，然后再逐組進行拼接。拼接后可利用屏蔽前的測序數(shù)據(jù)恢復contig中的重復序列，并利用正反向測序信息進行拼接和恢復正確性的檢驗。由于分組不會完全合理，還需要用blast等軟件對各contig進行比較以便去除冗余。然后進一步利用正反向測序信息恢復一些較長的重復序列，并構(gòu)建各contig之間的順序關(guān)系。這樣就基本完成了工作框架圖的拼接工作。
設(shè)計測序插入片段長度的分布由于在上述拼接過程中屏蔽掉了重復序列，這樣就會在基因組序列中留下一些洞(gap)。為了正確地填補這些洞，必須首先構(gòu)建正確的片段框架(scaffold)。在沒有詳細的物理圖譜等額外信息的情況下，就要精心設(shè)計測序插入片段的長度分布，為保證所有屏蔽重復序列后留下來的洞都能被適當長度的插入片段克隆覆蓋兩次以上。這就需要對插入片段長度分布進行設(shè)計。
插入片段覆蓋指定長度的洞的概率仍可用(3)式計算，只是為了構(gòu)建框架片段需要做一些小的修改。主要是要在洞的兩邊各留下50bp以上的序列，以便進行匹配識別。因此(3)式變?yōu)镻k=P(Yik=1)=CNk(L-F+100G)k(1-L-F+100G)N-k---(4)]]>上式表示起點為i，長度為F的洞被長度為L的插入片段覆蓋k次的概率。其中N為插入片段總數(shù)，G為基因組總長。若忽略基因組起始和結(jié)尾的區(qū)域，(4)式對任何點都成立。因此可略去下標i。
若要求洞被覆蓋兩次以上，概率為P2+＝1-P0-P1(5)由于當插入片段長度遠大于洞長時對構(gòu)建片段框架不利，而接近洞長時覆蓋效率又太低，根據(jù)經(jīng)驗，本發(fā)明規(guī)定長度L的片段只用于覆蓋0.2L至0.6L的洞。用人和水稻已有的序列進行檢驗，發(fā)現(xiàn)屏蔽重復序列后留下的最大洞在20-25kb左右，據(jù)此選定插入片段長度及其負責覆蓋洞長的數(shù)值見表2。
表2.插入片段長度及負責覆蓋洞的長度

測序成功率一般在90％左右，因此即使對所有插入片段都進行兩端測序，也仍會有20％左右克隆只能得到一端測序數(shù)據(jù)。所以在本發(fā)明的模型中規(guī)定20％測序片段沒有克隆另一端測序片段的數(shù)據(jù)。
如果沒有關(guān)于洞的長度分布的信息，本發(fā)明建議各長度插入片段的覆蓋度均取為10X，這樣對各種長度的洞覆蓋兩次的概率都在99％左右，就能取得較好的效果，見圖4。如果能根據(jù)已有的部分序列得到洞的長度分布，則可從上述初值出發(fā)，用公式(4)和(5)得到覆蓋各種長度洞的概率，再結(jié)合洞長分布可得到各長度洞未能被覆蓋的期望值，從而可計算出在這種插入片段分布下單位基因組長度上會留下多少洞。以此為目標函數(shù)，以總測序量為約束條件，可用非線性規(guī)劃的方法得到一定測序量下使遺留洞長達到最小的插入片段長度分布。一般情況下插入片段長度分布不需要這樣高的精度，也可采用excel等表格軟件輔助進行手動調(diào)整，其結(jié)果通常已可滿足使用需要。
權(quán)利要求
1.一種基于重復序列識別的全基因組測序數(shù)據(jù)的拼接方法，其特征在于該方法包括以下步驟(1)設(shè)定一個最小的DNA片段長度為15bp-20bp，計算散彈法測序數(shù)據(jù)中非重復片段出現(xiàn)的概率分布下列公式中各參數(shù)的含意為G基因組總長，L測序平均有效讀長，N成功測序反應數(shù)，F(xiàn)識別最小片段長度，定義一個隨機變量Yik描述用散彈法對全基因組測序中上述指定長度DNA片段出現(xiàn)K次的事件 若某點開始的片段出現(xiàn)次數(shù)為k，則有k個測序片段的起點應在基因組上[I-L+F，i]區(qū)間內(nèi)，而其它N-k個測序片段的起點不在此區(qū)間內(nèi)，這一區(qū)間的長度為L-F，若所有測序片段起點在基因組上隨機分布，則根據(jù)古典概型，上述隨機變量等于1的概率為P(Yik=1)=CNk(L-F+1G)k(1-L-F+1G)N-k---(1)]]>一次測序中出現(xiàn)次數(shù)為k的片段的平均個數(shù)可表示為E(Yk)=E(Σi=1GYik)=G·GNk(L-F+1G)k(1-L-F+1G)N-k----(2)]]>使用下式作為一次測序中出現(xiàn)次數(shù)為k的片段出現(xiàn)概率的估計值；Pk＝E(Yk)/G (3)(2)計算重復片段出現(xiàn)概率分布設(shè)片段為一個有m個拷貝的重復序列，出現(xiàn)在基因組中的m個不同位置，在散彈法測序數(shù)據(jù)集中的出現(xiàn)次數(shù)是所有位置出現(xiàn)次數(shù)的和，用Gmk表示有m個拷貝的重復序列在一次測序中的出現(xiàn)次數(shù)為k的概率，則上述關(guān)系可用數(shù)學公式表示為Gm0=P0m]]>Gm1=Cm1·P1·P0m-1]]>Gm2=Cm2·P12·P0m-2+Cm1·P2·P0m-1]]>Gm3=Cm3·P13·P0m-3+Cm2·C21·P1·P2·P0m-2+Cm1·P3·P0m-1]]>Gmj+＝1-Gm0-Gm1…-Gmj-1其中Gmj+表示出現(xiàn)次數(shù)為j和更多的概率；(3)重復序列的識別選取非重復片段出現(xiàn)概率最接近0.3％左右的次數(shù)為重復片段的判別標準，超過這一標準的片段就認為它屬于重復序列，否則就是非重復序列；(4)首先屏蔽重復序列，將上述散彈法測序數(shù)據(jù)中與識別出的片段相同的堿基改寫為N，屏蔽后剩余長度超過一個定值的測序數(shù)據(jù)仍進入拼接過程；(5)若目標基因組大小為1百萬-3千萬堿基，則屏蔽去重復序列后不分組直接進入拼接，若目標基因組大于上述范圍，則需要按照測序讀數(shù)之間的關(guān)連進行分組，然后進行拼接；(6)將得到的大片段中的N恢復成原有堿基，并利用同一克隆正反向測序的信息找出相關(guān)的大片段以及可能出現(xiàn)在它們之間的讀數(shù)，并將其連接起來；(7)所有能連接的片段都連接以后，再使用正反向測序信息把大片段排好順序，即得到目標基因組的工作框架圖。
全文摘要
本發(fā)明涉及一種基于重復序列識別的全基因組測序數(shù)據(jù)的拼接方法,首先計算散彈法測序數(shù)據(jù)中非重復片段和重復片段出現(xiàn)的概率分布,并根據(jù)這一概率分布確定重復序列的識別標準,然后用該標準屏蔽重復序列,再根據(jù)目標基因組的大小進行分組拼接,將得到的大片段中的N恢復成原有堿基,并利用同一克隆正反向測序的信息找出相關(guān)的大片段以及可能出現(xiàn)在它們之間的讀數(shù),并將其連接起來,所有能連接的片段都連接以后,再使用正反向測序信息把大片段排好順序,即得到目標基因組的工作框架圖。本發(fā)明的方法,具有提高效率、能處理復雜基因組、明顯減少錯誤拼接出現(xiàn)的概率、減少大量前期生物學實驗準備等優(yōu)點。
文檔編號C12P19/00GK1360057SQ0113485
公開日2002年7月24日 申請日期2001年11月16日 優(yōu)先權(quán)日2001年11月16日
發(fā)明者李松崗, 王俊, 蓋伊·王, 于軍, 汪建, 楊煥明, 倪培相, 韓玉軍, 黃顯剛, 張建國, 胡詠武 申請人:北京華大基因研究中心