專利名稱：短序列組裝中序列片段的過濾方法及系統(tǒng)的制作方法
技術(shù)領(lǐng)域：
本發(fā)明涉及基因工程技術(shù)領(lǐng)域，尤其涉及一種短序列組裝中序列片段的過濾方法及系統(tǒng)。
背景技術(shù)：
新測序技術(shù)產(chǎn)生的短序列有以下兩個特點第一，序列長度短；第二，數(shù)據(jù)量大。長序列組裝常用的Phrap等軟件均為基于序列間的交疊(overlap)來進行拼接組裝，此方法運用于短序列上會存在運算量太大的問題，沒有實際的應用價值。新興的短序列組裝受到內(nèi)存、時間等的限制，目前只在較小的原核生物基因組中成功應用。新一代測序分析存在以下難點第一，海量序列片段，基因組源序列的長度從十萬堿基(如豬痘病毒、大腸桿菌)到十億堿基(如黃種人、黃瓜、熊貓基因組)大小不等，而復雜環(huán)境(如海水、人體大腸等)宏基因組數(shù)據(jù)甚至會達到上百億堿基，而對這些樣本進行測序其覆蓋度需達到30倍到100倍，這使得產(chǎn)生的基因序列片段劇增，如亞洲黃種人的基因數(shù)據(jù)可達到ITB ;第二，短序列，隨著測序技術(shù)的發(fā)展，測序讀長呈不斷減小的趨勢，較第一代測序儀的測序長度顯著下降，例如454測序儀可以測到400bp, Sanger測序法的測序長度可達IOOObp到1200bp ;第三,測序錯誤，在測序產(chǎn)生序列片段的過程中可能伴隨由于熒光強度識別問題帶來測序誤差，例如有可能一個堿基T可能被測序儀讀出為A。這些錯誤是難以避免的，而且這個范圍通常是0. 5%到2%之間。這就意味著一個長度為75bp的源序列如果帶有1%的錯誤率，那么將導致有一半(1-(1-1%) 75=52. 9%)的測序產(chǎn)生序列片段可能有錯誤堿基。針對其中第二個問題，高通量的數(shù)據(jù)本身就可以生成大規(guī)模的k-mer節(jié)點，這些節(jié)點將被構(gòu)造成圖來分析，而由于測序錯誤的引入，將使得k-mer節(jié)點的數(shù)目增大5倍，例如人類基因組測序數(shù)據(jù)將會產(chǎn)生大約15G的k-mer ；由測序錯誤產(chǎn)生的k-mer，如果進入計算機進行直接處理，將會消耗巨大的內(nèi)存，例如人類基因組測序數(shù)據(jù)如果不進行序列過濾清洗的話，將會消耗大約2T的內(nèi)存來存儲這些k-mer所構(gòu)造的圖；測序數(shù)據(jù)中的錯誤序列還會在構(gòu)造的圖里面形成錯誤鏈接，Tip型錯誤，泡型錯誤，這些錯誤和源基因組序列本身的重復序列，基因突變點位等攪合在一起，這將使得后續(xù)的基因序列分析無法進行。因此，在短序列組裝前進行過濾，去除錯誤的k-mer，對序列的組裝和后續(xù)分析，尤其是大規(guī)模數(shù)據(jù)的分析，大基因組的組裝具有重要的意義。研究有效的序列過濾方法，節(jié)約內(nèi)存，提升計算性能成為一個亟待解決的問題。

發(fā)明內(nèi)容
本發(fā)明旨在解決上述現(xiàn)有技術(shù)中存在的問題，提出一種短序列組裝中序列片段的過濾方法，包括以下步驟接收測序序列；分別將接收到的測序序列逐個堿基滑動切割得到固定堿基長度的短串；將得到的所述短串的序列值及所述短串的頻率存儲為一個節(jié)點；計算所述短串頻率閾值；
將頻率小于閾值的短串過濾。優(yōu)選地,所述節(jié)點采用hash map存儲,其中，哈希鍵為所述序列值,值為所述節(jié)點。優(yōu)選地，所述將得到的所述短串的序列值及所述短串的頻率存儲為一個節(jié)點的步驟具體為根據(jù)當前節(jié)點的短串的序列值在已存儲的節(jié)點中查詢是否已存有當前節(jié)點；如果沒有查詢到當前節(jié)點，則添加所述當前節(jié)點；如果查詢到當前節(jié)點，則更新所述節(jié)點的頻率。

優(yōu)選地,所述節(jié)點中存儲短串和互補短串中序列值較大者或較小者。優(yōu)選地，所述閾值為T= 0 XCove, 0為分類模型參數(shù)，Covk為測序儀器設(shè)定的序列克隆倍數(shù)實際值。優(yōu)選地，所述計算所述短串頻率閾值中包括以下步驟以短串出現(xiàn)的頻率為橫坐標，以出現(xiàn)所述頻率的短串的個數(shù)為縱坐標，繪制頻率統(tǒng)計圖。優(yōu)選地，所述Covk的值為所述頻率統(tǒng)計圖上第一個波峰所在位置對應的覆蓋度。優(yōu)選地,所述Covk的計算方法步驟為a、對所有的短串按照出現(xiàn)頻率的個數(shù)排序,并把短串的個數(shù)按頻率的大小升序存入一個數(shù)組a中；b、刪除數(shù)組a中前面遞減的短串個數(shù)；C、用數(shù)組a的前j個數(shù)據(jù)求和來初始化SumO ；d、每次從數(shù)組a中取出第i個短串個數(shù)，加到Sumx里面，同時Sumx減去第i_j個頻率短串的個數(shù)，其中i大于j且i從j開始；e、如果Sunv^Sumx,回到步驟c,直到SunvPSunv進入下一步驟；f、用 j 除以 Sumx,即得到 Covk。本發(fā)明還提供了一種短序列組裝中序列片段的過濾系統(tǒng)，包括接收單元，用于接收測序序列；序列切割單元，用于分別將接收到的測序序列逐個堿基滑動切割得到固定堿基長度的短串；存儲統(tǒng)計單元，將得到的所述短串的序列值及所述短串的頻率存儲為一個節(jié)點；統(tǒng)計計算單元，用于計算所述短串頻率閾值；過濾單元，用于將頻率小于閾值的短串過濾。 優(yōu)選地,所述存儲統(tǒng)計單元包括查詢模塊，用于根據(jù)得到的短串的序列值在已存儲的節(jié)點中查詢是否已存有當前節(jié)點；節(jié)點添加模塊，用于在所述查詢模塊沒有查詢到當前節(jié)點時，添加當前節(jié)點；頻率更新模塊，用于在所述查詢模塊查詢到當前節(jié)點時，更新所述當前節(jié)點的頻率。本發(fā)明的有益效果在于，過濾了錯誤的短串，減小了組裝拼接的短串集合，減小了組裝拼接程序所需內(nèi)存，提高了組裝拼接程序的性能；在進行短串節(jié)點存儲的同時對短串出現(xiàn)的頻率進行了統(tǒng)計，操作簡單；誤差小。


圖1是本發(fā)明提供的序列片段的過濾方法的實現(xiàn)流程圖。圖2是本發(fā)明提供的序列片段的過濾的系統(tǒng)的結(jié)構(gòu)圖。圖3是本發(fā)明實施例中大腸桿菌的測序數(shù)據(jù)的短串頻率統(tǒng)計圖。圖4是本發(fā)明實施例中變異模型模擬測序數(shù)據(jù)的短串頻率統(tǒng)計圖。圖5是本發(fā)明實施例中454測序儀模型模擬測序數(shù)據(jù)的短串頻率統(tǒng)計圖。
具體實施例方式為了使本領(lǐng)域的技術(shù)人員更好的理解本申請的技術(shù)方案，下面將結(jié)合本申請實施例中的附圖，對本申請實施例中的技術(shù)方案進行清楚、完整的描述。應當理解，此處所描述的具體實施例僅僅用以解釋本發(fā)明，并不用于限定本發(fā)明。在本發(fā)明的實施例中，通過分別將接收到的測序序列逐個堿基滑動切割得到固定堿基長度的短串(k-mer)，并將得到的各短串的序列值存儲，統(tǒng)計得到的各所述短串出現(xiàn)的頻率，繪制所述短串的頻率統(tǒng)計圖，計算所述短串頻率閾值，將頻率小于閾值的短串過濾。圖1所示為本發(fā)明實施例提供的短序列組裝中序列片段過濾方法的實現(xiàn)流程，詳述如下在步驟SlOl中，接收測序序列；在步驟S102中，分別將接收到的測序序列逐個堿基滑動切割得到固定堿基長度的短串(k-mer)；在步驟S103中，將得到的所述短串的序列值及所述短串的頻率存儲為一個節(jié)點；在步驟S104中，計算所述短串頻率閾值；在步驟S105中，將頻率小于閾值的短串過濾。在本發(fā)明的實施例中，測序序列的堿基長度為25-75，切割成固定堿基長度為21-31的短串。然而，切割得到的短串的長度小于測序序列的長度，其長度可以根據(jù)測序序列的長度和實際情況設(shè)定。每個節(jié)點存儲相應短串的序列值和頻率。這里，可采用longIongint類型文件存儲所述節(jié)點,其存儲格式如下[seq :64, frequency :16,...]；其中，seq存儲短串的序列值，所述序列值的計算方法是使用2位存儲一個核苷酸序列，如A用00表不,G用01表不,C用10表不,T用11表不,順序編碼下去生成一個占64位的整數(shù)值，并且，考慮到對于偶數(shù)長度的短串，其互補短串可能為它本身，例如短串GATC的互補短串為GATC本身。為了防止這種混淆，短串的長度均為奇數(shù)，另外，由于本發(fā)明實施例中數(shù)據(jù)結(jié)構(gòu)的限制，短串的長度取21-31的奇數(shù)frequency用16位存儲所述短串出現(xiàn)的次數(shù)，即頻率，頻率的取值范圍為
;其后面的位數(shù)還可以用來存儲其他值，例如，可以存儲刪除標記closed，以標識所述短串是否被刪除；也可以存儲使用標記in_use，以標識所述短串是否被使用過，還可以存儲其他標識。上述步驟S103具體為步驟1，根據(jù)當前節(jié)點的短串的序列值在已存儲的節(jié)點中查詢是否已存有當前節(jié)占.步驟2，如果沒有查詢到當前節(jié)點，則添加所述當前節(jié)點；
步驟3，如果查詢到當前節(jié)點，則更新所述當前節(jié)點的頻率。本發(fā)明在存儲各節(jié)點的同時，對短串的頻率進行了統(tǒng)計。在本發(fā)明的實施例中，使用hash map存儲各節(jié)點，哈希鍵為序列值，值為節(jié)點。例如序列為AAAAA的短串(其互補序列為TTTTT)，其序列值為1111111111，頻率初始值為1，將其序列值1111111111作為鍵在hash map中查詢是否已經(jīng)存有當前節(jié)點,如果沒有查詢到當前節(jié)點,則添加所述當前節(jié)點存儲到hash map中，其值為所述短串的序列值1111111111，頻率初始值為I ;如果查詢到當前節(jié)點，則對所述當前節(jié)點頻率進行更新，增加I。完成后，執(zhí)行步驟2，查找下一個短串，直至完成全部短串的查找。為了降低存儲節(jié)點所需的空間，作為本發(fā)明的一個優(yōu)選實施例，只用一個節(jié)點存儲互補的兩個短串，節(jié)點的序列值取互補的兩個短串中較大的序列值。如果一個短串的序列值小于其互補短串的序列值，則節(jié)點存儲所述互補短串的序列值，例如上例中序列AAAAA的序列值存的就是其互補短串TTTTT的值；如果一個短串的序列值大于其互補短串的序列值，則節(jié)點存儲所述短串的序列值。當然，節(jié)點的序列值也可以存儲互補的兩個短串中較小的序列值。當然，也可以用其他結(jié)構(gòu)對各節(jié)點進行存儲，例如可以用樹結(jié)構(gòu)進行存儲，使用hash map存儲各節(jié)點在內(nèi)存和使用上與用樹狀結(jié)構(gòu)存儲近似,但是用hash map存儲各節(jié)點在訪問和修改速度上都明顯優(yōu)于樹結(jié)構(gòu)。步驟S104計算所述短串頻率閾值，在本實施例中頻率閾值的計算方法如下所述閾值為T= 0 XCove, 0為分類模型參數(shù)，Covk為測序儀器設(shè)定的序列克隆倍數(shù)的實際值。分類模型參數(shù)的范圍一般在0-10%，當分類模型參數(shù)偏小時，被過濾的短串(k-mer)較少,可能保留了更多的錯誤k_mer ;當分類模型參數(shù)偏大時,被過濾的短串(k-mer)較多,可能會勿將正確的k_mer也過濾掉了，對后續(xù)序列拼接組裝或基因分析造成影響。因此，分類模型參數(shù)根據(jù)實際計算的內(nèi)存條件，后續(xù)序列拼接所使用算法特點等因素進行選擇。測序儀器設(shè)定的序列克隆倍數(shù)是一個理論值，在實際測序過程中可以設(shè)定為某一固定值，但是，由于測序儀的誤差和測序過程中的操作誤差，測序儀器設(shè)定的序列克隆倍數(shù)的實際值與理論值相差較大，因此，要根據(jù)測序結(jié)果對其重新進行計算。在本發(fā)明的一個實施例中，以短串出現(xiàn)的頻率為橫坐標，出現(xiàn)所述頻率的短串的個數(shù)為縱坐標繪制頻率統(tǒng)計圖。根據(jù)上述的頻率統(tǒng)計圖，所述Covk的值為所述頻率統(tǒng)計圖上第一個波峰所在位置對應的覆蓋度。例如，選取大腸桿菌的測序數(shù)據(jù)進行k-mer頻率統(tǒng)計，所述頻率統(tǒng)計圖橫坐標為短串出現(xiàn)的頻率，縱坐標為出現(xiàn)所述頻率的短串的個數(shù)，結(jié)果如圖3所示，第一個波峰所對應的點為(62，12. 68)，從圖3可讀出Covk值為62。在本發(fā)明的另一個實施例中，所述Covk的值可按如下步驟進行計算a、對所有的短串按照出現(xiàn)頻率的個數(shù)排序，并把短串的個數(shù)按頻率的大小升序存入一個數(shù)組a中；b、刪除數(shù)組a中前面遞減的短串個數(shù)；C、用數(shù)組a的前j個數(shù)據(jù)求和來初始化SumO ；d、每次從數(shù)組a中取出第i個短串個數(shù)，加到Sumx里面，同時Sumx減去第i_j個頻率短串的個數(shù)，其中i大于j且i從j開始；e、如果Sunv^Sumx,回到步驟c,直到SunvASumx,進入下一步驟；f、用 j 除以 Sumx,即得到 Covk。通過設(shè)定的分類模型參數(shù)和計算出的測序儀器設(shè)定的序列克隆倍數(shù)實際值，可以得到頻率閾值，將頻率小于閾值的短串過濾。本領(lǐng)域的普通技術(shù)人員可以理解，實現(xiàn)上述實施例方法中的全部或部分步驟是可以通過程序來指令相關(guān)的硬件來完成，所述的程序可以在存儲于一計算機可讀取存儲介質(zhì)中，所述的存儲介質(zhì)可以為R0M/RAM、磁盤、光盤等，所述程序用來執(zhí)行以下步驟I，接收測序序列；2，分別將接收到的測序序列逐個堿基滑動切割得到固定堿基長度的短串(k-mer)； 3，將得到的所述短串的序列值及所述短串的頻率存儲為一個節(jié)點；4，計算所述短串頻率閾值；5，將頻率小于閾值的短串過濾。圖2所示為本發(fā)明實施例提供的短序列組裝中序列片段過濾的系統(tǒng)的結(jié)構(gòu)，為了便于說明僅示出了與本發(fā)明實施例相關(guān)的部分。所述短序列組裝中序列片段過濾的系統(tǒng)可以用于短序列組裝或基因分析中，其中接收單元201，用于接收測序序列。序列切割單元202，用于分別將接收到的測序序列逐個堿基滑動切割得到固定堿基長度的短串，其實現(xiàn)方式如上所述，在此不再一一贅述。存儲統(tǒng)計單元203，用于將得到的所述短串的序列值及所述短串的頻率存儲為一個節(jié)點，其實現(xiàn)方式如上所述，在此不再一一贅述。統(tǒng)計計算單元204，用于計算所述短串頻率閾值。過濾單元205，用于將頻率小于閾值的短串過濾。其中，所述存儲統(tǒng)計單元203包括查詢模塊2031，用于根據(jù)得到的短串的序列值在已存儲的節(jié)點中查詢是否已存有當前節(jié)點。節(jié)點添加模塊2032，用于在所述查詢模塊沒有查詢到當前節(jié)點時，添加當前節(jié)點，其實現(xiàn)方式如上所述，不再一一贅述。頻率統(tǒng)計模塊2033，用于在所述查詢模塊查詢到當前節(jié)點時，更新所述節(jié)點的頻率，所述節(jié)點頻率增加I。以下結(jié)合具體的測序儀器模擬數(shù)據(jù)對本發(fā)明的過濾系統(tǒng)進行誤差分析。首先利用變異模型生成的模擬測序數(shù)據(jù)進行驗證。變異模型假設(shè)一個短序列中每個位置測序儀出錯的可能性相同。令RefSeq的長度為N，并且RefSeq中重疊(r印eats)所占的比例為￠,測序儀器的誤差設(shè)定為a，k為de novo拼接算法中所設(shè)定的k_mer的長度。于是，理論上可以得到正確k-mer的個數(shù)為Kptjsitive，錯誤k-mer的個數(shù)為Knegative,計算公式分別為
權(quán)利要求
1.一種短序列組裝中序列片段的過濾方法，其特征在于，所述方法包括以下步驟 接收測序序列； 分別將接收到的測序序列逐個堿基滑動切割得到固定堿基長度的短串； 將得到的所述短串的序列值及所述短串的頻率存儲為一個節(jié)點； 計算所述短串頻率閾值； 將頻率小于閾值的短串過濾。
2.根據(jù)權(quán)利要求1所述的過濾方法，其特征在于，所述節(jié)點采用hashmap存儲,其中，哈希鍵為所述序列值，值為所述節(jié)點。
3.根據(jù)權(quán)利要求1所述的過濾方法，其特征在于，所述將得到的所述短串的序列值及所述短串的頻率存儲為一個節(jié)點的步驟具體為 根據(jù)當前節(jié)點的短串的序列值在已存儲的節(jié)點中查詢是否已存有當前節(jié)點； 如果沒有查詢到當前節(jié)點，則添加所述當前節(jié)點； 如果查詢到當前節(jié)點，則更新所述當前節(jié)點的頻率。
4.根據(jù)權(quán)利要求1所述的過濾方法，其特征在于，所述節(jié)點中存儲短串和互補短串中序列值較大者或較小者。
5.根據(jù)權(quán)利要求1所述的過濾方法，其特征在于，所述閾值為T=0 XCove, 0為分類模型參數(shù)，Cove為測序儀器設(shè)定的序列克隆倍數(shù)實際值。
6.根據(jù)權(quán)利要求5所述的過濾方法，其特征在于，所述計算所述短串頻率閾值的步驟包括以下步驟以短串出現(xiàn)的頻率為橫坐標，以出現(xiàn)所述頻率的短串的個數(shù)為縱坐標，繪制頻率統(tǒng)計圖。
7.根據(jù)權(quán)利要求6所述的過濾方法，其特征在于，所述Covk的值為所述頻率統(tǒng)計圖上第一個波峰所在位置對應的覆蓋度。
8.根據(jù)權(quán)利要求5所述的過濾方法，其特征在于，所述Covk的計算方法步驟為 a、對所有的短串按照出現(xiàn)頻率的個數(shù)排序，并把短串的個數(shù)按出現(xiàn)頻率的大小升序存入一個數(shù)組a中； b、刪除數(shù)組a中前面遞減的短串個數(shù)； C、用數(shù)組a的前j個數(shù)據(jù)求和來初始化SumO ； d、每次從數(shù)組a中取出第i個短串個數(shù)，加到Sumx里面，同時Sumx減去第i_j個頻率短串的個數(shù)，其中i大于j且i從j開始； e、如果Sunv^Sumx,回到步驟c,直到SunvASunv進入下一步驟； f、用j除以Sumx,即得到Covk。
9.一種短序列組裝中序列片段的過濾系統(tǒng)，其特征在于，所述系統(tǒng)包括 接收單元，用于接收測序序列； 序列切割單元，用于分別將接收到的測序序列逐個堿基滑動切割得到固定堿基長度的短串； 存儲統(tǒng)計單元，將得到的所述短串的序列值及所述短串的頻率存儲為一個節(jié)點； 統(tǒng)計計算單元，用于計算所述短串頻率閾值； 過濾單元，用于將頻率小于閾值的短串過濾。
10.根據(jù)權(quán)利要求9所述的系統(tǒng)，其特征在于，所述存儲統(tǒng)計單元包括查詢模塊，用于根據(jù)得到的短串的序列值在已存儲的節(jié)點中查詢是否已存有當前節(jié)點節(jié)點添加模塊，用于在所述查詢模塊沒有查詢到當前節(jié)點時，添加當前節(jié)點；頻率更新模塊，用于在所述查詢模塊查詢到當前節(jié)點時，更新所述當前節(jié)點的頻率。
全文摘要
本發(fā)明公開了一種短序列組裝中序列片段的過濾方法，包括以下步驟接收測序序列；分別將接收到的測序序列逐個堿基滑動切割得到固定堿基長度的短串；將得到的所述短串的序列值及所述短串的出現(xiàn)頻率存儲為一個節(jié)點；計算所述短串頻率閾值；將頻率小于閾值的短串過濾。本發(fā)明還提供了短序列組裝中序列片段的過濾系統(tǒng)。本發(fā)明的有益效果在于，過濾了錯誤的短串，減小了組裝拼接的短串集合，減小了組裝拼接程序所需內(nèi)存，提高了組裝拼接程序的性能；在進行短串節(jié)點存儲的同時對短串出現(xiàn)的頻率進行了統(tǒng)計，操作簡單；誤差小。
文檔編號G06F19/22GK103065067SQ20121057572
公開日2013年4月24日 申請日期2012年12月26日 優(yōu)先權(quán)日2012年12月26日
發(fā)明者孟金濤, 魏彥杰, 曾理, 成杰峰, 馮圣中 申請人:深圳先進技術(shù)研究院