本發(fā)明涉及生物信息技術(shù)領(lǐng)域，具體來(lái)說(shuō)，涉及一種基因組序列的組裝方法、結(jié)構(gòu)變異探測(cè)方法和相應(yīng)的系統(tǒng)。

背景技術(shù)：

基因組測(cè)序是開展分子生物學(xué)研究的重要技術(shù)。通過(guò)對(duì)一個(gè)物種的基因組進(jìn)行測(cè)序，研究人員可以獲得這個(gè)物種的基因組堿基序列，它作為這個(gè)物種的遺傳序列模版，為基因、轉(zhuǎn)錄、調(diào)控、修飾等層面進(jìn)行定性或定量的研究，探索生命現(xiàn)象背后的分子機(jī)制提供了重要參照。完成測(cè)序后，通過(guò)將被測(cè)物種的基因組與其他物種的基因組進(jìn)行比較，研究人員可以發(fā)現(xiàn)它們?cè)诨蚪M水平上的差異，這為揭示遺傳變異、自然或人工選擇的機(jī)制提供了信息，從而為優(yōu)質(zhì)基因的篩選、物種的改良培育提供了指導(dǎo)。此外，基因組測(cè)序還可以幫助尋找多倍體物種的雜合位點(diǎn)或雜合區(qū)段，是研究雜合性與生命現(xiàn)象的關(guān)系的重要基礎(chǔ)。

第二代測(cè)序技術(shù)是目前應(yīng)用的最廣泛的測(cè)序技術(shù)。和第一代測(cè)序技術(shù)相比，它具有通量高、成本低的特點(diǎn)。在第二代測(cè)序技術(shù)發(fā)展的最初階段，所測(cè)到的序列長(zhǎng)度比較短(堿基對(duì)數(shù)目通常為30-40bp)而且堿基辨識(shí)的質(zhì)量不夠高；隨著技術(shù)的不斷改進(jìn)，目前能夠測(cè)量的序列長(zhǎng)度大幅度增加(超過(guò)100bp)，同時(shí)堿基辨識(shí)質(zhì)量也有了很大的改進(jìn)。第二代測(cè)序的一個(gè)重要的特點(diǎn)是，它可以從兩端對(duì)一個(gè)很長(zhǎng)的片段進(jìn)行測(cè)序，得到這個(gè)長(zhǎng)片段兩端的堿基序列，因此使用第二代測(cè)序技術(shù)可以獲得高通量的雙末端測(cè)序序列。

將測(cè)序序列組裝成基因組是計(jì)算生物領(lǐng)域的基本問(wèn)題。因?yàn)闇y(cè)序儀所能測(cè)量的序列長(zhǎng)度遠(yuǎn)小于基因組長(zhǎng)度，所以在測(cè)序后需要對(duì)所有測(cè)序序列進(jìn)行組裝，推斷它們的相對(duì)位置，還原出被測(cè)的基因組。組裝基因組面臨著以下幾個(gè)挑戰(zhàn)：(1)第二代測(cè)序技術(shù)的數(shù)據(jù)具有很高的通量，大量的數(shù)據(jù)會(huì)增加組裝的時(shí)間和計(jì)算設(shè)備上被占用的存儲(chǔ)空間；(2)基因組上有很多相似度很高、或是重復(fù)出現(xiàn)的區(qū)段，它們的存在給推測(cè)測(cè)序序列的相對(duì)位置增加了很大的不確定性；(3)對(duì)于雜合度高的基因組，需要組裝出不同的倍型，同時(shí)還要確定不同倍型的位置關(guān)系，找到雜合區(qū)域。

現(xiàn)有的基因組組裝方法在原理上主要分為兩類。一類是基于debruijn圖的方法，該方法的主要操作是：對(duì)于每一個(gè)測(cè)序序列，每隔一個(gè)堿基切割出一個(gè)特定長(zhǎng)度的子序列(通常稱為k-mer，k表示子序列的長(zhǎng)度)；利用所有被切割出來(lái)的子序列構(gòu)造debruijn圖；進(jìn)行一定的糾錯(cuò)操作后，在圖上尋找路徑，每條路徑被推斷為被測(cè)基因組上的片段。這種方法對(duì)雜合度低，重復(fù)度低的基因組組裝效果會(huì)比較好；而對(duì)于雜合度高，重復(fù)度高的基因組不是很理想。此外，基于debruijn圖的方法不容易給出測(cè)序序列之間的疊落關(guān)系，不利于進(jìn)行統(tǒng)計(jì)評(píng)估；同時(shí)，所切割出的子序列長(zhǎng)度明顯短于測(cè)序序列，會(huì)降低特異性，導(dǎo)致在圖上尋找路徑時(shí)出現(xiàn)錯(cuò)誤。另一類方法是基于測(cè)序序列疊落關(guān)系的，該方法對(duì)每?jī)蓷l序列進(jìn)行比對(duì)，根據(jù)全部比對(duì)結(jié)果推斷序列的疊落關(guān)系。這類方法是以測(cè)序序列為單位的，而不是k-mer，容易從序列水平進(jìn)行統(tǒng)計(jì)評(píng)估。然而，這一類方法需要對(duì)每?jī)蓷l測(cè)序序列進(jìn)行比對(duì)，時(shí)間復(fù)雜度高，對(duì)測(cè)序通量小的第一代測(cè)序技術(shù)可以應(yīng)用，但不適用于高通量的第二代測(cè)序技術(shù)。同時(shí)，對(duì)于重復(fù)度高的基因組，這類方法的效果也不一定理想，會(huì)出現(xiàn)拷貝數(shù)減少；對(duì)于雜合度比較高的區(qū)域，當(dāng)基因組上同一位置的兩個(gè)倍型差異比較大時(shí)，來(lái)自兩個(gè)倍型的測(cè)序序列不容易被整合到一個(gè)疊陣中，導(dǎo)致倍型的丟失。

針對(duì)相關(guān)技術(shù)中的上述問(wèn)題，目前尚未提出有效的解決方案。

技術(shù)實(shí)現(xiàn)要素：

針對(duì)相關(guān)技術(shù)中的上述問(wèn)題，本發(fā)明提出一種基因組序列的組裝方法、變異探測(cè)方法和組裝系統(tǒng)，能夠?qū)崿F(xiàn)了序序列的高效連接，實(shí)現(xiàn)測(cè)序序列的基因組裝。

本發(fā)明的技術(shù)方案是這樣實(shí)現(xiàn)的：

根據(jù)本發(fā)明的一個(gè)方面，提供了一種基因組序列的組裝方法。

該組裝方法包括：

(1)通過(guò)預(yù)定的映射算法將樣品的被測(cè)基因組的測(cè)序序列映射到參考基因組，得到單映射測(cè)序序列疊陣集，其中，樣品的測(cè)序序列為利用高通量測(cè)序技術(shù)測(cè)得，參考基因組已知并與樣品的基因組相近；

(2)基于經(jīng)過(guò)預(yù)處理的參考基因組對(duì)單映射測(cè)序序列疊陣集中的測(cè)序序列進(jìn)行篩選，所得篩選結(jié)果根據(jù)覆蓋度再次篩選，得到篩選后的單映射測(cè)序序列疊陣集；

(3)通過(guò)單方向測(cè)序序列信息對(duì)篩選后的單映射測(cè)序序列疊陣集進(jìn)行切割，得到初始預(yù)組裝疊陣集，將當(dāng)前預(yù)組裝疊陣集的初始值設(shè)置為初始預(yù)組裝疊陣集；

(4)確定當(dāng)前預(yù)組裝疊陣集中每個(gè)疊陣的相對(duì)位置，形成組裝基因組架構(gòu)；

(5)對(duì)組裝基因組架構(gòu)中的每個(gè)疊陣進(jìn)行延拓，得到每個(gè)疊陣的一致序列；

(6)將組裝基因組架構(gòu)中的相鄰疊陣的一致序列中符合預(yù)定連接規(guī)則的一致序列進(jìn)行連接，得到樣品的當(dāng)前的組裝基因組；

(7)根據(jù)被測(cè)基因組上同源序列的差異調(diào)整預(yù)定的映射算法的映射參數(shù)，通過(guò)調(diào)整后的該預(yù)定的映射算法將樣品的被測(cè)基因組的測(cè)序序列映射到當(dāng)前的組裝基因組，得到當(dāng)前預(yù)組裝疊陣集；

對(duì)當(dāng)前預(yù)組裝疊陣集迭代執(zhí)行步驟(4)、(5)和(6)，迭代次數(shù)為任何非負(fù)整數(shù)。

其中，步驟(1)中的預(yù)定的映射算法中包括預(yù)定的映射參數(shù)，預(yù)定的映射參數(shù)包括以下至少之一：

被測(cè)基因組與參考基因組的差異預(yù)期；

被測(cè)基因組的長(zhǎng)度、測(cè)序序列的長(zhǎng)度和質(zhì)量特征；

其中，映射參數(shù)用于提供判別任意一個(gè)測(cè)序序列和參考基因組上某個(gè)位置起始的子序列的相似度是否達(dá)到測(cè)序序列成功映射到參考基因組位置的準(zhǔn)則。

此外，該步驟(1)包括：

在將樣品的被測(cè)基因組的測(cè)序序列映射到參考基因組后，將被測(cè)基因組中映射到參考基因組上多個(gè)位置的測(cè)序序列去除，得到單映射測(cè)序序列疊陣集。

另外，該步驟(2)中對(duì)參考基因組的預(yù)處理包括：

對(duì)參考基因組進(jìn)行自映射，得到參考基因組中的若干唯一性序列區(qū)域。

此外，在執(zhí)行步驟(1)對(duì)測(cè)序序列進(jìn)行映射時(shí)的映射率低于預(yù)定標(biāo)準(zhǔn)的情況下，則進(jìn)行下述操作：

在執(zhí)行步驟(2)后，對(duì)于篩選后的測(cè)序序列疊陣集，在每一個(gè)位置，選擇最大頻數(shù)的堿基，用最大頻數(shù)的堿基更新參考基因組的唯一性序列區(qū)域上對(duì)應(yīng)位置的堿基；

調(diào)整預(yù)定的映射算法的預(yù)定的映射參數(shù)，基于已經(jīng)更新過(guò)唯一性序列區(qū)域的參考基因組，重新執(zhí)行步驟(1)和步驟(2)。

其中，在執(zhí)行步驟(1)的映射操作和步驟(2)的篩選操作時(shí)，如果被測(cè)基因組的雙末端測(cè)序序列數(shù)據(jù)集的一對(duì)同伴序列的兩端都被單映射到當(dāng)前預(yù)組裝疊陣集，則一對(duì)同伴序列的映射坐標(biāo)信息在步驟(4)中用于形成組裝基因組架構(gòu)；

如果被測(cè)基因組的雙末端測(cè)序序列數(shù)據(jù)集的一對(duì)同伴序列中的至少一端被單映射到當(dāng)前預(yù)組裝疊陣集，則一對(duì)同伴序列中的至少一端的映射坐標(biāo)信息在步驟(5)中用于疊陣延拓；

其中，雙末端測(cè)序序列數(shù)據(jù)集包括多個(gè)具有不同庫(kù)長(zhǎng)的片段庫(kù)。

此外，該步驟(3)中的對(duì)篩選后的單映射測(cè)序序列疊陣集進(jìn)行切割包括：

對(duì)于參考基因組上的每一個(gè)堿基，計(jì)算覆蓋該堿基的所有左向測(cè)序序列的尾長(zhǎng)的最大值w1，以及計(jì)算覆蓋該堿基的所有右向測(cè)序序列的尾長(zhǎng)的最大值w2；

如果w1或者w2小于一個(gè)預(yù)定的閾值w，則將該堿基標(biāo)記為切割位點(diǎn)。預(yù)定的閾值w為整數(shù)，且0≤w≤lmax，lmax為所有測(cè)序序列長(zhǎng)度的最大值；

其中，左向測(cè)序序列的尾長(zhǎng)和右向測(cè)序序列的尾長(zhǎng)的定義包括：

對(duì)于參考基因組上的任一個(gè)堿基，該堿基將覆蓋該堿基的每一條測(cè)序序列分成左右兩部分；其中，如果左側(cè)部分的長(zhǎng)度大于右側(cè)部分的長(zhǎng)度，則稱該測(cè)序序列為左向測(cè)序序列，并且右側(cè)部分的長(zhǎng)度為該左向測(cè)序序列的尾長(zhǎng)；如果右側(cè)部分的長(zhǎng)度大于左側(cè)部分的長(zhǎng)度，則稱該測(cè)序序列為右向測(cè)序序列，并且左側(cè)部分的長(zhǎng)度為該右向測(cè)序序列的尾長(zhǎng)。

另外，該步驟(4)包括：

利用測(cè)序序列的庫(kù)長(zhǎng)信息以及單映的測(cè)序序列在當(dāng)前預(yù)組裝疊陣集中的坐標(biāo)，確定當(dāng)前預(yù)組裝疊陣集中任意兩個(gè)疊陣間的距離范圍；

對(duì)當(dāng)前預(yù)組裝疊陣集中的所有疊陣進(jìn)行排列，使每?jī)蓚€(gè)疊陣間的距離與確定的對(duì)應(yīng)該每?jī)蓚€(gè)疊陣間的距離范圍相匹配。

此外，該步驟(5)包括：

在當(dāng)前預(yù)組裝疊陣集中的每個(gè)疊陣的每個(gè)端點(diǎn)附近設(shè)定一個(gè)范圍，利用單映的測(cè)序序列在疊陣中的坐標(biāo)信息，確定疊陣?yán)锓秶鷥?nèi)的測(cè)序序列的同伴序列，同伴序列與疊陣的一致序列共同構(gòu)成從端點(diǎn)向外延拓疊陣的測(cè)序信息庫(kù)，預(yù)定范圍與疊陣中的測(cè)序序列所屬的片段庫(kù)的庫(kù)長(zhǎng)一致；

對(duì)測(cè)序信息庫(kù)中的所有序列按照局部比對(duì)的算法進(jìn)行比對(duì)，得到兩兩對(duì)比結(jié)果；

利用圖論的深度優(yōu)先算法整合兩兩比對(duì)結(jié)果，形成每個(gè)端點(diǎn)附近延拓后的疊陣；

基于延拓后的疊陣定義延拓的一致序列。

此外，該步驟(6)包括：

利用局部比對(duì)算法判斷相鄰疊陣一致序列是否存在重疊情況；

在存在重疊的情況下，將該相鄰的疊陣的一致序列進(jìn)行連接，得到樣品的當(dāng)前的組裝基因組。

此外，該步驟(7)中的預(yù)定的映射算法中包括預(yù)定的映射參數(shù)，預(yù)定的映射參數(shù)包括以下至少之一：

被測(cè)基因組上同源序列的差異；

被測(cè)基因組的長(zhǎng)度、測(cè)序序列的長(zhǎng)度和質(zhì)量特征；

映射參數(shù)用于提供判別任意一個(gè)測(cè)序序列和當(dāng)前基因組上某個(gè)位置起始的子序列的相似度是否達(dá)到測(cè)序序列成功映射到參考基因組位置的準(zhǔn)則。

該步驟(7)包括：

在將樣品的被測(cè)基因組的測(cè)序序列映射到當(dāng)前的組裝基因組后，將被測(cè)基因組中映射到當(dāng)前的組裝基因組上多個(gè)位置的測(cè)序序列去除，得到當(dāng)前單映射測(cè)序序列疊陣集；

對(duì)當(dāng)前單映射測(cè)序序列疊陣集執(zhí)行步驟(3)，得到當(dāng)前預(yù)組裝疊陣集。

此外，基于上述組裝方法中的任意一種組裝方法還可進(jìn)行雙倍體序列的組裝。

根據(jù)本發(fā)明的另一方面，提供了一種應(yīng)用上述組裝方法中任意一種組裝方法的結(jié)構(gòu)變異探測(cè)方法。

該結(jié)構(gòu)變異探測(cè)方法包括：

根據(jù)組裝方法對(duì)不同樣本的基因組之間的結(jié)構(gòu)變異情況進(jìn)行探測(cè)，探測(cè)的信息包括步驟(3)中對(duì)單映射測(cè)序序列疊陣集進(jìn)行切割時(shí)形成的斷點(diǎn)信息。

根據(jù)本發(fā)明的再一方面，提供了一種基因組序列的組裝系統(tǒng)。

該組裝系統(tǒng)包括：

映射模塊，用于通過(guò)預(yù)定的映射算法將樣品的被測(cè)基因組的測(cè)序序列映射到參考基因組，得到單映射測(cè)序序列疊陣集，其中，樣品的測(cè)序序列為利用高通量測(cè)序技術(shù)測(cè)得，參考基因組已知并與樣品的基因組相近；

篩選模塊，用于基于經(jīng)過(guò)預(yù)處理的參考基因組對(duì)單映射測(cè)序序列疊陣集中的測(cè)序序列進(jìn)行篩選，所得篩選結(jié)果根據(jù)覆蓋度再次篩選，得到篩選后的單映射測(cè)序序列疊陣集；

切割模塊，用于通過(guò)單方向測(cè)序序列信息對(duì)篩選后的單映射測(cè)序序列疊陣集進(jìn)行切割，得到初始預(yù)組裝疊陣集；

架構(gòu)模塊，用于確定初始預(yù)組裝疊陣集中每個(gè)疊陣的相對(duì)位置，形成組裝基因組架構(gòu)；

延拓模塊，用于對(duì)組裝基因組架構(gòu)中的每個(gè)疊陣進(jìn)行延拓，得到每個(gè)疊陣的一致序列；

連接模塊，用于將組裝基因組架構(gòu)中的相鄰疊陣的一致序列中符合預(yù)定連接規(guī)則的一致序列進(jìn)行連接，得到樣品的當(dāng)前的組裝基因組；

調(diào)整映射模塊，用于根據(jù)被測(cè)基因組上同源序列的差異調(diào)整預(yù)定的映射算法的預(yù)定的映射參數(shù)，通過(guò)調(diào)整的該預(yù)定的映射算法將樣品的被測(cè)基因組的測(cè)序序列映射到當(dāng)前的組裝基因組，得到當(dāng)前預(yù)組裝疊陣集；

延拓模塊和連接模塊進(jìn)一步用于對(duì)調(diào)整映射模塊中的當(dāng)前預(yù)組裝疊陣集進(jìn)行操作。

本發(fā)明通過(guò)將被測(cè)基因組的測(cè)序序列和參考基因組進(jìn)行映射，并對(duì)映射結(jié)果進(jìn)行切割，以及將切割后的疊陣進(jìn)行組裝和延拓，從而實(shí)現(xiàn)了測(cè)序序列的高效連接，實(shí)現(xiàn)測(cè)序序列的基因組裝。

附圖說(shuō)明

為了更清楚地說(shuō)明本發(fā)明實(shí)施例或現(xiàn)有技術(shù)中的技術(shù)方案，下面將對(duì)實(shí)施例中所需要使用的附圖作簡(jiǎn)單地介紹，顯而易見地，下面描述中的附圖僅僅是本發(fā)明的一些實(shí)施例，對(duì)于本領(lǐng)域普通技術(shù)人員來(lái)講，在不付出創(chuàng)造性勞動(dòng)的前提下，還可以根據(jù)這些附圖獲得其他的附圖。

圖1是根據(jù)本發(fā)明實(shí)施例的基因組序列的組裝方法的總體流程圖；

圖2是根據(jù)本發(fā)明實(shí)施例的基因組組裝方法的詳細(xì)流程圖；

圖3是根據(jù)本發(fā)明實(shí)施例的確定參考基因組上自映射唯一性區(qū)域的示意圖；

圖4是根據(jù)本發(fā)明實(shí)施例的根據(jù)唯一性條件進(jìn)行篩選的示意圖；

圖5是根據(jù)本發(fā)明實(shí)施例的利用單方向測(cè)序序列信息確定疊陣切割位點(diǎn)方法的示意圖；

圖6是根據(jù)本發(fā)明實(shí)施例的利用雙末端信息連接疊陣估計(jì)相鄰疊陣距離的示意圖；

圖7是根據(jù)本發(fā)明實(shí)施例的收集用于延拓疊陣一端所需要的測(cè)序序列的方法示意圖；

圖8是根據(jù)本發(fā)明實(shí)施例的計(jì)算兩個(gè)測(cè)序序列的最優(yōu)位移的示意圖；

圖9是根據(jù)本發(fā)明實(shí)施例的對(duì)疊陣進(jìn)行延拓的方法的示意圖；

圖10(a)～圖10(d)是根據(jù)本發(fā)明實(shí)施例的連接相鄰疊陣的一致序列示意圖；

圖11是根據(jù)本發(fā)明實(shí)施例的通過(guò)再映射和局部組裝獲取雙倍體序列方法示意圖。

圖12是根據(jù)本發(fā)明實(shí)施例的基因組序列的組裝系統(tǒng)的框圖。

具體實(shí)施方式

下面將結(jié)合本發(fā)明實(shí)施例中的附圖，對(duì)本發(fā)明實(shí)施例中的技術(shù)方案進(jìn)行清楚、完整地描述，顯然，所描述的實(shí)施例僅僅是本發(fā)明一部分實(shí)施例，而不是全部的實(shí)施例?；诒景l(fā)明中的實(shí)施例，本領(lǐng)域普通技術(shù)人員所獲得的所有其他實(shí)施例，都屬于本發(fā)明保護(hù)的范圍。

根據(jù)本發(fā)明的實(shí)施例，提供了一種基因組序列的組裝方法。

為了便于理解本發(fā)明實(shí)施例的基因組序列的組裝方法，下面先對(duì)本發(fā)明實(shí)施例的相關(guān)技術(shù)名詞作出如下定義和解釋：

基因組：包含一個(gè)生物體所有遺傳信息的遺傳序列，由四種堿基排列而成，通常表示為含有a，c，g，t四種字符的序列，每一種字符代表一種堿基。

雙倍體：大多生物個(gè)體體細(xì)胞含有兩套基因組，一套來(lái)自于父本，另一套來(lái)自于母本；這兩套基因組基本一致，但也存在差異；存在差異的位置為雜合位點(diǎn)或者雜合區(qū)域；在雜合位點(diǎn)或者雜合區(qū)域，兩套基因組序列不一致，具有兩個(gè)倍型。

被測(cè)基因組：需要確定堿基序列的基因組。

測(cè)序序列：由測(cè)序平臺(tái)輸出的某一生物體基因組子序列的測(cè)量結(jié)果，由若干個(gè)代表四種堿基的字符(a，c，g，t)組成，代表基因組的一個(gè)子序列。測(cè)序序列的一端為5’端，另一端為3’端；通常在測(cè)序平臺(tái)給出的文件中，測(cè)序序列的左端為5’端，右端為3’端。

bp：英語(yǔ)basepair的縮寫。100bp指一個(gè)堿基序列有100個(gè)堿基。

第一代測(cè)序技術(shù)：最早出現(xiàn)的測(cè)序技術(shù)，利用若干生物化學(xué)技術(shù)，將被測(cè)基因組序列切割成多個(gè)子序列，并確定每個(gè)子序列堿基組成；測(cè)序序列長(zhǎng)度可以到500bp左右，但產(chǎn)生的測(cè)序序列數(shù)據(jù)量較小。

第二代測(cè)序技術(shù)：也稱作新一代測(cè)序技術(shù)，英文名稱為nextgenerationsequencing，英文縮寫為ngs；其特點(diǎn)是能并行地產(chǎn)生大量的雙末端測(cè)序序列數(shù)據(jù)，每一條測(cè)序序列的長(zhǎng)度大多為100bp-200bp。

雙末端測(cè)序序列：第二代測(cè)序技術(shù)在進(jìn)行操作時(shí)，將被測(cè)的序列切割成若干片段，片段的長(zhǎng)度可小可大；之后在片段的兩端進(jìn)行測(cè)序，在一端得到一個(gè)測(cè)序序列，在另一端得到一個(gè)測(cè)序序列；如果片段過(guò)長(zhǎng)，中間部分一般無(wú)法測(cè)到；所得到的兩個(gè)測(cè)序序列對(duì)應(yīng)于同一個(gè)片段，稱為雙末端測(cè)序序列。

片段庫(kù)、庫(kù)長(zhǎng)：在測(cè)序平臺(tái)給出的同一個(gè)雙末端測(cè)序序列文件中，所有雙末端測(cè)序序列所屬的片段長(zhǎng)度相近，測(cè)序平臺(tái)會(huì)估計(jì)出一個(gè)平均值；在下文中，稱測(cè)序平臺(tái)給出的每一個(gè)雙末端測(cè)序序列文件為片段庫(kù)，稱其估計(jì)出片段長(zhǎng)度平均值為庫(kù)長(zhǎng)。

堿基誤讀：由于測(cè)序技術(shù)的誤差使得測(cè)序序列上某些堿基不同于真實(shí)堿基。

質(zhì)量值：在測(cè)序過(guò)程中，測(cè)序序列的每一個(gè)堿基都有被誤讀的可能性，質(zhì)量值是這一可能性的反應(yīng)。測(cè)序序列的每一個(gè)堿基都對(duì)應(yīng)一個(gè)這樣的質(zhì)量值，質(zhì)量值越高，被誤讀的可能性越低。

參考基因組：一個(gè)已完成測(cè)序或者組裝，明確了堿基序列的基因組，在研究過(guò)程中被用作模板或者參照。

插入：指相對(duì)于參考基因組，測(cè)序序列的某兩個(gè)相鄰的堿基之間額外插入了一段堿基序列。

刪失：指相對(duì)于參考基因組，測(cè)序序列丟失一段或若干段堿基序列。

indel：一個(gè)插入或者刪失。

測(cè)序序列映射：對(duì)于一個(gè)測(cè)序序列，在參考基因組上尋找一個(gè)與其長(zhǎng)度大體一致的子序列，該子序列與測(cè)序序列完全一樣，或者存在差異，但差異在預(yù)先設(shè)定的標(biāo)準(zhǔn)之內(nèi)。通常以堿基替換個(gè)數(shù)，indel長(zhǎng)度來(lái)衡量二者的差異。

堿基替換：測(cè)序序列被映射到參考基因組上后某些堿基和參考基因組上對(duì)應(yīng)的堿基不相同。

映射結(jié)果：對(duì)于一個(gè)測(cè)序序列，如果參考基因組上存在子序列，且子序列與測(cè)序序列的差異在設(shè)定的標(biāo)準(zhǔn)之內(nèi)，則該測(cè)序序列映射成功。映射結(jié)果包括：(1)參考基因組上的子序列最左端堿基的坐標(biāo)，該值作為測(cè)序序列在參考基因組上的映射坐標(biāo)；(2)測(cè)序序列的映射方向，可以按5’至3’的方向映射，也可以按3’至5’的方向映射；(3)測(cè)序序列與參考基因組上子序列的比對(duì)信息，包括兩者之間存在的堿基替換，indel情況。當(dāng)參考基因組上存在多個(gè)與測(cè)序序列差異在設(shè)定標(biāo)準(zhǔn)之內(nèi)的子序列時(shí)，該測(cè)序序列具有多個(gè)映射結(jié)果，每一個(gè)映射結(jié)果對(duì)應(yīng)于參考基因組上的一個(gè)子序列。

基因組組裝：利用測(cè)序序列還原被測(cè)基因組，組裝后得到的基因組稱為組裝基因組。

測(cè)序序列疊陣(可以簡(jiǎn)稱為疊陣)：英文術(shù)語(yǔ)為contig或layout；表示由若干測(cè)序序列形成的集合，在該集合當(dāng)中，每?jī)蓷l測(cè)序序列之間的疊落關(guān)系或者相對(duì)距離是確定的。測(cè)序序列疊陣可以用含有a，c，g，t，ф五種字符的矩陣表示；矩陣的某個(gè)元素為ф表示矩陣的該位置上不存在堿基；矩陣的每一列對(duì)應(yīng)于疊陣的每一個(gè)位點(diǎn)且每一列至少含有一個(gè)非ф字符；矩陣的列數(shù)為疊陣的長(zhǎng)度；矩陣的每一行代表一個(gè)測(cè)序序列，每行第一個(gè)非ф字符所在列的列標(biāo)為該行所代表的測(cè)序序列在疊陣中的坐標(biāo)。

疊陣集：由疊陣組成的集合。

架構(gòu)(或組裝基因組架構(gòu))：英文術(shù)語(yǔ)為scaffold；它由一系列按照一定順序排列后的疊陣組成，其中所有疊陣的相對(duì)位置關(guān)系已知，作為組裝基因組的一個(gè)框架。

一致序列：英文術(shù)語(yǔ)為consensus；對(duì)于疊陣的每一列，根據(jù)其包含的所有堿基信息推斷出一個(gè)堿基，并將所推斷出的堿基連接后得到的序列；推斷堿基的一個(gè)可行的方法為，對(duì)疊陣的每一列，在其包含的所有堿基中取頻率最大者。

在清楚了相關(guān)技術(shù)名詞的定義和解釋后，在本發(fā)明的一個(gè)實(shí)施例中，為使得基因組組裝能夠高效地進(jìn)行，同時(shí)提高具有高雜合度，或者高重復(fù)度的基因組的組裝準(zhǔn)確性，提出了一種基于高通量測(cè)序技術(shù)的基因組組裝方法。

具體的，本發(fā)明的實(shí)施例中提出的基因組組裝方法以雙末端測(cè)序序列和一個(gè)參考基因組作為輸入。輸入的參考基因組可以是來(lái)自相近物種的已被測(cè)序的基因組；也可以是使用其他組裝方法對(duì)測(cè)序序列進(jìn)行組裝，但效果不夠理想、需要進(jìn)行修正的組裝基因組。其中，在一個(gè)實(shí)施例中，雙末端測(cè)序序列來(lái)自被測(cè)基因組，雙末端測(cè)序序列可以包含多個(gè)具有不同庫(kù)長(zhǎng)的片段庫(kù)。鑒于被測(cè)基因組與參考基因組來(lái)自相同或者相近物種，二者具有一定的相似度，本發(fā)明借助參考基因組，利用測(cè)序序列映射工具，實(shí)現(xiàn)被測(cè)基因組的組裝。如圖1所示，本發(fā)明包括下述步驟：

步驟1、將測(cè)序序列映射至參考基因組，并經(jīng)過(guò)后續(xù)操作得到初始預(yù)組裝疊陣集；所述初始預(yù)組裝疊陣集為一個(gè)疊陣集，它是步驟2進(jìn)行組裝的基礎(chǔ)；

其中，在一個(gè)實(shí)施例中，可根據(jù)參考基因組、測(cè)序序列的特征，測(cè)序序列所屬基因組與參考基因組差異大小的先驗(yàn)信息，以及對(duì)映射的靈敏度、特異度的預(yù)期，設(shè)計(jì)映射所用的參數(shù)；

其中，測(cè)序序列的特征指測(cè)序序列的長(zhǎng)度分布情況和質(zhì)量值分布情況；二者所屬物種基因組差異大小指兩個(gè)基因組之間出現(xiàn)snp、indel的頻率以及indel的長(zhǎng)度分布情況；而對(duì)于庫(kù)長(zhǎng)不同，或者測(cè)序序列長(zhǎng)度不同的片段庫(kù)，可以設(shè)計(jì)不同的映射參數(shù)；

在一個(gè)實(shí)施例中，在完成映射后，將映射到參考基因組上多個(gè)位置的測(cè)序序列去除，獲得單映射測(cè)序序列疊陣集；

另外，在一個(gè)實(shí)施例中，對(duì)參考基因組進(jìn)行自映射，獲得參考基因組上具有自映射唯一性的區(qū)域，并利用唯一性條件對(duì)所述單映射測(cè)序序列疊陣進(jìn)行篩選，得到篩選后的疊陣集；如果映射率不夠理想，可以根據(jù)映射率的具體情況調(diào)整映射參數(shù)重新進(jìn)行映射，也可以基于篩選后得到的疊陣集對(duì)參考基因組進(jìn)行更新，并重新執(zhí)行步驟1；獲得篩選的疊陣集之后，評(píng)估其中每個(gè)疊陣的連續(xù)性，在必要的位置進(jìn)行切割，得到初始預(yù)組裝疊陣集；

步驟2、基于步驟1所得初始預(yù)組裝疊陣集進(jìn)行基因組的組裝，得到組裝基因組；

利用雙末端序列的庫(kù)長(zhǎng)信息估計(jì)初始預(yù)組裝疊陣集中每?jī)蓚€(gè)疊陣的距離，重新排列疊陣使每?jī)蓚€(gè)疊陣的距離與估計(jì)值吻合，得到組裝基因組架構(gòu)；對(duì)于每個(gè)疊陣，將位于其兩端的測(cè)序序列的同伴序列，與位于疊陣兩端的測(cè)序序列進(jìn)行比對(duì)，使疊陣向兩端延拓，并推斷一致序列；比對(duì)相鄰的疊陣的一致序列，根據(jù)比對(duì)結(jié)果精確連接一致序列，得到組裝基因組。

步驟3、將測(cè)序序列映射至組裝基因組，獲得當(dāng)前預(yù)組裝疊陣集，基于當(dāng)前預(yù)組裝疊陣集重新執(zhí)行步驟2，實(shí)現(xiàn)迭代組裝。步驟3可以執(zhí)行，也可以不執(zhí)行。

完成上述步驟1，步驟2和步驟3后，將組裝基因組輸出至標(biāo)準(zhǔn)格式的文件。

以下詳細(xì)闡述本發(fā)明的方法和原理。圖2詳細(xì)地示出了本發(fā)明實(shí)施例的相關(guān)處理流程。

本發(fā)明實(shí)施例的步驟1根據(jù)參考基因組，雙末端測(cè)序序列的具體情況，以及對(duì)二者差異的先驗(yàn)知識(shí)設(shè)計(jì)映射參數(shù)，將測(cè)序序列映射至參考基因組，獲得相對(duì)于參考基因組的單映的測(cè)序序列疊陣集，并利用唯一性條件對(duì)疊陣進(jìn)行篩選，得到篩選后的疊陣集，再對(duì)其中的疊陣進(jìn)行必要的切割，得到初始預(yù)組裝疊陣集。其具體包括以下步驟：

步驟11、設(shè)計(jì)映射參數(shù)。映射所用的主要為：

a:測(cè)序序列與參考基因組之間的錯(cuò)誤匹配個(gè)數(shù)上界m。該參數(shù)由參考基因組和測(cè)序序列所屬物種的差異率γ，測(cè)序序列的長(zhǎng)度l，以及測(cè)序序列堿基誤讀的平均比率。一般地，可以將m設(shè)為測(cè)序序列長(zhǎng)度的5％或6％。如果預(yù)測(cè)參考基因組和測(cè)序序列所屬物種的差異較大，可以將m調(diào)整為測(cè)序序列長(zhǎng)度的10％至15％。

b：可以被探測(cè)的indel的長(zhǎng)度的最大值max_indel，該參數(shù)反映參考基因組和測(cè)序序列所屬物種的差異；同時(shí)，該參數(shù)也會(huì)影響映射的速度。一般地，max_indel可以選擇為5。

c：參數(shù)s，表示一個(gè)測(cè)序序列最多可以被映射到參考基因組上s個(gè)位置。如果一個(gè)測(cè)序序列被映射到了參考基因組上超過(guò)s個(gè)位置，即認(rèn)為該測(cè)序序列映射失敗。優(yōu)選地，可以將s設(shè)為大于10的整數(shù)。

如果采用單個(gè)子序列-延拓的映射算法，還需設(shè)計(jì)以下參數(shù)：

d：完全匹配子序列的長(zhǎng)度下限值k。該參數(shù)表示在映射時(shí)需要在測(cè)序序列上找到一個(gè)長(zhǎng)度不小于k的子序列，同時(shí)在參考基因組上能夠找到一個(gè)與之完全相同的子序列，以該子序列在參考基因組上的位置作為測(cè)序序列的初步定位。該參數(shù)的選取依賴于參考基因組和測(cè)序序列所屬物種的差異率，測(cè)序序列的長(zhǎng)度，測(cè)序序列堿基辨識(shí)的平均錯(cuò)誤率，以及參考基因組的總長(zhǎng)度。

e：每個(gè)測(cè)序序列的子序列搜尋數(shù)目上限值u。該參數(shù)設(shè)得越大，所尋找的子序列數(shù)目越多，但也會(huì)增加計(jì)算時(shí)間。優(yōu)選地，可以將其設(shè)為20。

在設(shè)計(jì)參數(shù)時(shí)，可以對(duì)同一個(gè)片段庫(kù)的所有測(cè)序序列設(shè)計(jì)共同的參數(shù)。此時(shí)，可以以所有測(cè)序序列長(zhǎng)度的眾數(shù)l代替l。對(duì)于測(cè)序序列長(zhǎng)度分布，或者堿基辨識(shí)質(zhì)量值分布具有明顯差異的不同片段庫(kù)，可以設(shè)置不同的參數(shù)。

步驟12、對(duì)于每一個(gè)片段庫(kù)，利用步驟11所設(shè)計(jì)的參數(shù)，將該片段庫(kù)的雙末端測(cè)序序列映射至參考基因組，將映射到多個(gè)位置的測(cè)序序列從映射結(jié)果中去除，得到單映射測(cè)序序列疊陣集。如果映射率低于預(yù)定的標(biāo)準(zhǔn)，則進(jìn)入步驟13，否則結(jié)束步驟1的操作。

步驟13、利用唯一性條件對(duì)單映測(cè)序序列疊陣進(jìn)行篩選，得到篩選后的疊陣集。其具體包括以下步驟：

步驟131、對(duì)參考基因組進(jìn)行自映射，以獲得參考基因組上具有自映射唯一性的區(qū)段。對(duì)于長(zhǎng)度分布具有明顯差異的不同序列片段庫(kù)，可以分別進(jìn)行操作，得到不同自映射唯一性區(qū)域。步驟131具體包括以下步驟：

步驟1311、從參考基因組上每隔一個(gè)堿基截取一個(gè)長(zhǎng)度為l的序列，將這些序列映射到參考基因組上。優(yōu)選地，可以利用步驟11中設(shè)定的參數(shù)完成此步的映射。

步驟1312、對(duì)參考基因組上的每一個(gè)堿基，設(shè)定一個(gè)深度值，將所有堿基的深度值設(shè)定為0。遍歷步驟1311中所有序列的映射結(jié)果，進(jìn)行以下操作：如果一條序列映射到了參考基因組上唯一的位置，則將參考基因組上該序列所覆蓋住的所有堿基的深度值加1；如果該序列被映射到了參考基因組上多于一個(gè)位置，則對(duì)于每一個(gè)成功映射的位置，該序列所覆蓋的所有堿基的深度值加1。

步驟1313、記錄參考基因組上所有具有自映射唯一性的區(qū)域。所述自映射唯一性的區(qū)域?yàn)閰⒖蓟蚪M上的一個(gè)區(qū)間，滿足以下條件：該區(qū)間中所有堿基的深度值均等于l；任意包含該區(qū)間的其它區(qū)間均含有深度值不等于l的堿基。

圖3為確定參考基因組上自映射唯一性區(qū)域的方法示意圖，如圖3所示，完成自映射后，所有連續(xù)的且深度為l的堿基形成一個(gè)自映射唯一性區(qū)域。

步驟132、遍歷步驟12中得到的所有單映的測(cè)序序列，按照和步驟1312相同的方法重新計(jì)算參考基因組上每一個(gè)堿基的深度值；取dα為深度值分布的上α分位數(shù)，其中α為大于0小于0.5的數(shù)，優(yōu)選地，以將其取為0.05；

步驟133、對(duì)于步驟12中得到的每一個(gè)單映的測(cè)序序列，檢驗(yàn)其是否同時(shí)滿足以下兩個(gè)唯一性條件。如果其不能同時(shí)滿足，則將其從所處的單映射測(cè)序序列疊陣中去除，最后得到一系列新的疊陣，稱這些疊陣組成的疊陣集為“經(jīng)過(guò)唯一性條件篩選后的單映射測(cè)序序列疊陣集”；

其中，第一唯一性條件：在參考基因組上，所述測(cè)序序列完全包含于具有自映射唯一性的區(qū)域；第二唯一性條件：步驟132中計(jì)算出的參考基因組上被所述測(cè)序序列覆蓋的所有堿基的深度值小于dα。

圖4為根據(jù)唯一性條件進(jìn)行篩選的示意圖，如圖4所示，不滿足唯一性條件的測(cè)序序列從單映測(cè)序序列疊陣中被去除。

步驟14、如果步驟12的映射率、或者單映射的測(cè)序序列所占比率低于預(yù)定的標(biāo)準(zhǔn)，可以進(jìn)行以下操作；兩個(gè)操作可以全部進(jìn)行，也可以擇一進(jìn)行：

操作一：調(diào)整映射參數(shù)；優(yōu)選地，可以將測(cè)序序列與參考基因組之間的錯(cuò)誤匹配個(gè)數(shù)上界m調(diào)大；返回步驟12；

操作二：設(shè)定頻率值下界θ和深度值上界d，對(duì)于參考基因組上的每個(gè)堿基，基于步驟133得到的疊陣，計(jì)算該堿基的深度值以及覆蓋住該堿基位點(diǎn)的a，c，g，t四種核苷酸的頻率；如果該堿基的深度值大于d，且最大的頻率值qmax大于θ，則將參考基因組上該堿基替換為qmax所對(duì)應(yīng)的核苷酸，返回步驟11。

步驟15、利用單方向測(cè)序序列信息評(píng)估步驟133所得疊陣集中的疊陣的連續(xù)性，在必要位置對(duì)疊陣進(jìn)行切割，得到初始預(yù)組裝疊陣集。

在具體闡述步驟15的操作之前，引入下述定義：

任取參考基因組上的一個(gè)堿基，考查覆蓋該堿基的每一條測(cè)序序列，堿基將該測(cè)序序列分成左右兩部分。如果左側(cè)部分的長(zhǎng)度大于右側(cè)，則稱該測(cè)序序列為“左向測(cè)序序列”，同時(shí)稱右側(cè)部分的長(zhǎng)度為該左向測(cè)序序列的尾長(zhǎng)；如果右側(cè)部分的長(zhǎng)度大于左側(cè)，則稱該測(cè)序序列為“右向測(cè)序序列”，同時(shí)稱左側(cè)部分的長(zhǎng)度為該右向測(cè)序序列的尾長(zhǎng)。

步驟15具體包括以下步驟：

步驟151、對(duì)于參考基因組上的每一個(gè)堿基，計(jì)算覆蓋該堿基的所有左向測(cè)序序列的尾長(zhǎng)的最大值，記之為w1；同時(shí)計(jì)算覆蓋該堿基的所有右向測(cè)序序列的尾長(zhǎng)的最大值，記之為w2.如果w1或者w2小于一個(gè)預(yù)定的閾值w，則將所述堿基標(biāo)記為切割位點(diǎn)。所述w為一整數(shù)，0≤w≤lmax，lmax為所有測(cè)序序列長(zhǎng)度的最大值。

圖5以參考基因組上的堿基b為例示出了w1和w2的計(jì)算方式。圖中的星號(hào)代表每條測(cè)序序列的中點(diǎn)，如果該中點(diǎn)位于b的左側(cè)，則該測(cè)序序列為一左向測(cè)序序列，反之則為一右向測(cè)序序列。

步驟152、在步驟151所得的所有切割位點(diǎn)處，對(duì)步驟133得到的疊陣集中的疊陣進(jìn)行切割，得到切割后的疊陣集。具體的操作為：在參考基因組上從左至右掃描所有公共切割位點(diǎn)，對(duì)每一個(gè)切割位點(diǎn)，從步驟133得到的疊陣集中找到覆蓋該切割位點(diǎn)的疊陣；將該疊陣分割為兩個(gè)疊陣，其中一個(gè)包含該疊陣中所有映射至所述切割位點(diǎn)左側(cè)的測(cè)序序列，另一個(gè)包含該疊陣中所有映射至所述切割位點(diǎn)右側(cè)的測(cè)序序列；分割之后得到的每個(gè)疊陣中，測(cè)序序列的相對(duì)位置關(guān)系與其在進(jìn)行分割之前的疊陣中的相對(duì)位置關(guān)系一致。

下文中，至步驟3之前，若無(wú)特殊說(shuō)明，出現(xiàn)的所有“疊陣”二字均指步驟15所得疊陣集，即切割之后得到初始預(yù)組裝疊陣集中的疊陣。

本發(fā)明的步驟2基于步驟1得到的初始預(yù)組裝疊陣集實(shí)現(xiàn)基因組組裝，得到組裝基因組。其具體包括以下步驟：

步驟21、利用單映測(cè)序序列在疊陣上的坐標(biāo)信息以及雙末端測(cè)序序列的庫(kù)長(zhǎng)信息估計(jì)任意兩個(gè)疊陣的距離，并將所有疊陣進(jìn)行排列，使每?jī)蓚€(gè)疊陣的距離與估計(jì)值相吻合。優(yōu)選地，兩個(gè)疊陣x1和x2的距離可以采用步驟211和步驟212進(jìn)行估計(jì)：

步驟211、對(duì)兩個(gè)疊陣中的測(cè)序序列進(jìn)行掃描，如果x1中的一個(gè)測(cè)序序列r1與x2中的一個(gè)測(cè)序序列r2為雙末端測(cè)序序列，則計(jì)算r1的左端到x1的右端的距離，以及r2的右端到x2的左端的距離；用r1和r2所在的片段庫(kù)的庫(kù)長(zhǎng)減去兩個(gè)距離之和，得到一個(gè)差值，該差值作為x1與x2的距離的一個(gè)觀測(cè)值。

步驟212、統(tǒng)計(jì)步驟211得到的x1與x2的距離的觀測(cè)值個(gè)數(shù)，如果個(gè)數(shù)大于設(shè)定的下界，則取這些觀測(cè)值的中位數(shù)作為x1與x2的距離的估計(jì)值。

完成每?jī)蓚€(gè)疊陣的距離估計(jì)之后，將所有疊陣進(jìn)行排列，使得排列后的每?jī)蓚€(gè)疊陣的距離與估計(jì)值相吻合。

圖6為步驟21所述方法的示意圖。如圖6所示，i，ii和iii為三個(gè)疊陣，每條虛線連接的兩個(gè)箭頭表示一對(duì)雙末端測(cè)序序列。根據(jù)雙末端信息，將疊陣i，ii和iii排列為i→iii→ii，并估計(jì)相鄰疊陣之間的距離，得到組裝基因組架構(gòu)。

步驟22、采用基于測(cè)序序列兩兩比對(duì)的方法在左右兩端延拓每個(gè)疊陣。首先闡述將一個(gè)疊陣向右延拓的具體步驟，向左延拓可以用類似的方法實(shí)現(xiàn)。將一個(gè)疊陣x向右延拓的具體操作為：

步驟221、基于單映射測(cè)序序列在疊陣上的坐標(biāo)信息，收集將向右延拓疊陣x時(shí)所需的測(cè)序序列，計(jì)算每個(gè)測(cè)序序列對(duì)于x的先驗(yàn)坐標(biāo)，具體包括以下步驟：

步驟2211、建立集合set并將其初始化為空集。記lenx為疊陣x的長(zhǎng)度。以x的最左端為坐標(biāo)原點(diǎn)，計(jì)算x中所有測(cè)序序列的起始坐標(biāo)。

步驟2212、遍歷x中全體測(cè)序序列，對(duì)每一個(gè)測(cè)序序列進(jìn)行以下操作：

記r為所述測(cè)序序列，記posr為其在x上的起始坐標(biāo)，記insert_sizer為r所屬片段庫(kù)的庫(kù)長(zhǎng)。檢驗(yàn)r是否滿足以下兩個(gè)條件：(1)lenx-t·insert_sizer≤posr，t為預(yù)先設(shè)定的不小于1的數(shù)；(2)r在x上從左至右為5’端至3’端。如果r同時(shí)滿足上述兩個(gè)條件，則將r的同伴序列r′的反向互補(bǔ)序列添加到set中，以posr+insert_sizer-lr′作為r′對(duì)于x的先驗(yàn)坐標(biāo)，其中l(wèi)r′為r′的長(zhǎng)度。。

步驟2213、將位于x右端的一部分序列添加到set中，以這些序列在x上的坐標(biāo)作為先驗(yàn)坐標(biāo)。

所得集合set保存所有向右延拓疊陣x所需的測(cè)序序列。

圖7為步驟221的一個(gè)簡(jiǎn)單示意圖。圖中所有實(shí)線箭頭表示靠近疊陣x右端的測(cè)序序列，其中方向向右的測(cè)序序列的同伴序列表示為虛線箭頭，這些同伴序列即為所要收集用于延拓x的測(cè)序序列。

步驟222、設(shè)定兩個(gè)測(cè)序序列重疊部分的堿基匹配數(shù)目下界match_bound(優(yōu)選地將其設(shè)定為大于20的整數(shù))；設(shè)定所述匹配數(shù)目與重疊部分長(zhǎng)度的比值下界ratio_bound(優(yōu)選地將其設(shè)定為大于0.9且小于1的數(shù))；設(shè)定先驗(yàn)坐標(biāo)之差的閾值pos_bound(該閾值為不大于測(cè)序序列長(zhǎng)度的非負(fù)數(shù))。對(duì)于set中任意兩個(gè)測(cè)序序列r1和r2，如果二者先驗(yàn)坐標(biāo)之差的絕對(duì)值小于pos_bound，則計(jì)算r1相對(duì)于r2的最優(yōu)位移，以及r2相對(duì)于r1的相對(duì)位移。

記r1的長(zhǎng)度為記r2的長(zhǎng)度為所述r1相對(duì)于r2的最優(yōu)位移按照以下步驟計(jì)算：

步驟2221、設(shè)定三個(gè)變量s，max_ratio和optimal_shift(r1，r2)。將s初始化為將max_ratio初始化為ratio_bound，將optimal_shift初始化為正無(wú)窮。

步驟2222、計(jì)算兩個(gè)數(shù)值overlap(r1，r2，s)和match(r1，r2，s)，具體的計(jì)算方法如下：

當(dāng)s＜0時(shí)，

當(dāng)s≥0時(shí)，

和的定義方式如下：如果r1[i-s]＝r2[i]，則否則如果r1[i]＝r2[i+s],則否則

步驟2223、如果match(r1,r2,s)≥match_bound,則計(jì)算比值

ratio＝match(r1，r2，s)/overlap(r1，r2，s)。

如果ratio≥max_ratio，則將max_ratio替換為ratio，將optimal_shift替換為s。

步驟2224、如果則結(jié)束操作，否則用s+1代替s，返回步驟2222。

完成上述操作后得到的optimal_shift(r1，r2)即為測(cè)序序列r1相對(duì)于r2的最優(yōu)位移，其可以為正無(wú)窮。

圖8為計(jì)算最優(yōu)位移的一個(gè)實(shí)施例。如圖8所示，若將r2相對(duì)于r1向右平移5bp，二者的不匹配數(shù)目為15；若將r2相對(duì)于r1向右平移10bp，二者的不匹配數(shù)目為9；若將r2相對(duì)于r1向右平移8bp，二者的不匹配數(shù)目為0。因此r2相對(duì)于r1的最優(yōu)位移為8，r1相對(duì)于r2的最優(yōu)位移為-8。

步驟223、構(gòu)建有向圖g，g的每一個(gè)結(jié)點(diǎn)為集合set中每一個(gè)測(cè)序序列；對(duì)于r中任意兩條測(cè)序序列r1和r2，如果optimal_shift(r1，r2)為負(fù)，則在g中添加一條由r1指向r2的邊；如果optimal_shift(r1，r2)為正且不等于正無(wú)窮，則在g中添加一條由r2指向r1的邊。

步驟224、在步驟223構(gòu)建的有向圖g中，找出入度為0且先驗(yàn)坐標(biāo)最小的結(jié)點(diǎn)。如果g中不存在入度為0的結(jié)點(diǎn)，則找出先驗(yàn)坐標(biāo)最小的測(cè)序序列所對(duì)應(yīng)的結(jié)點(diǎn)。以所述被找出的結(jié)點(diǎn)作為初始結(jié)點(diǎn)，進(jìn)行以下操作：

步驟2241、對(duì)g中所有結(jié)點(diǎn)設(shè)定訪問(wèn)狀態(tài)，將初始結(jié)點(diǎn)的訪問(wèn)狀態(tài)設(shè)定為已訪問(wèn)，將其余結(jié)點(diǎn)的訪問(wèn)狀態(tài)設(shè)定為未訪問(wèn)；

步驟2242、從初始結(jié)點(diǎn)起，按照深度優(yōu)先原則對(duì)g進(jìn)行遍歷；在遍歷過(guò)程中，對(duì)于每一個(gè)被訪問(wèn)到的結(jié)點(diǎn)，將其訪問(wèn)狀態(tài)更改為已訪問(wèn)，并在該結(jié)點(diǎn)所指向的所有未被訪問(wèn)的結(jié)點(diǎn)中，選取相對(duì)于該結(jié)點(diǎn)最優(yōu)位移最小的一個(gè)進(jìn)行下一步的訪問(wèn)；如果該結(jié)點(diǎn)的出度為0，或者該結(jié)點(diǎn)所指向的所有結(jié)點(diǎn)均已被訪問(wèn)，則根據(jù)遍歷過(guò)程記錄從初始結(jié)點(diǎn)到該結(jié)點(diǎn)的路徑；

步驟2243、在步驟2242記錄的所有路徑中，選出含結(jié)點(diǎn)數(shù)最多的一條。

步驟225、利用步驟2243選出的路徑中所包含的測(cè)序序列構(gòu)造一個(gè)疊陣，記這個(gè)疊陣為y，具體的構(gòu)造方式如下：

記所述路徑所包含的結(jié)點(diǎn)依次為r1，…，rn，其中n為結(jié)點(diǎn)總數(shù)。將r1在y中的坐標(biāo)設(shè)為1；對(duì)于任意的正整數(shù)i(2≤i≤n)，將ri-1在y中的坐標(biāo)與optimal_shift(ri，ri-1)相加，作為ri在y中的坐標(biāo)。

步驟226、整合x和y，將x向右延拓，推斷一致序列。其具體操作包括以下步驟：

步驟2261、如果y中含有來(lái)自x的測(cè)序序列，則尋找一個(gè)同時(shí)屬于x和y的測(cè)序序列，記它在x和y中的坐標(biāo)分別為coorx和coory；對(duì)于y中的每一個(gè)測(cè)序序列，將其在y中的坐標(biāo)與(coorx-coory)相加，作為其在延拓后的疊陣x中的坐標(biāo)；對(duì)于同時(shí)包含于集合set以及x中的測(cè)序序列，如果其不在y中，則將其從x中去掉；如果y中不包含來(lái)自x的測(cè)序序列，則以x自身作為延拓后的疊陣。

步驟2262、推斷x的一致序列，記作c(x)。

上述步驟221至步驟226為向右延拓疊陣x，推斷得到一條或兩條一致序列的方法。向左延拓x時(shí)，可以將x整體取反向互補(bǔ)，之后同樣進(jìn)行步驟221至步驟226的操作，最后將所得到的一條或者兩條一致序列取反向互補(bǔ)。

在具體實(shí)施過(guò)程中，可以并行地對(duì)每一個(gè)疊陣執(zhí)行步驟22的操作。

步驟23、連接相鄰兩個(gè)疊陣的一致序列，得到組裝基因組。

首先闡述兩條序列seq1和seq2，且seq1在左，seq2在右的連接方法。一個(gè)優(yōu)選的方案如下：

截取seq1不同長(zhǎng)度的后綴序列，以及seq2不同長(zhǎng)度的前綴序列。如果seq1的某個(gè)后綴序列suffix_seq1，和seq2的某個(gè)前綴序列prefix_seq2能夠完全匹配，或者僅存在很少的堿基替換或者拆入/刪失，則認(rèn)為seq1和seq2能夠成功連接，否則認(rèn)為二者不能成功連接。如果seq1和seq2能夠成功連接，則將seq2中位于prefix_seq2右側(cè)的部分接到seq1右端，形成連接后的序列，記prefix_seq2的長(zhǎng)度為二者公共部分的長(zhǎng)度；如果seq1和seq2不能成功連接，則在seq1右端接上若干字符n，將seq2接在這些n的右端，得到將seq1和seq2連接后的序列。seq1和seq2連接后的序列記作seq1οseq2。

從左至右順次掃描組裝基因組架構(gòu)，對(duì)每?jī)蓚€(gè)相鄰、且已被延拓的疊陣的一致序列按照上述方案進(jìn)行連接，得到組裝基因組。

步驟3、如果步驟2所得組裝基因組不夠理想，可以以其作為參考基因組，將測(cè)序序列進(jìn)行映射，獲得當(dāng)前預(yù)組裝疊陣集，并基于當(dāng)前預(yù)組裝疊陣集執(zhí)行步驟2，實(shí)現(xiàn)迭代組裝，從而提高參考基因組的評(píng)估指標(biāo)。其具體包括以下步驟：

步驟31、根據(jù)同源序列的差異和測(cè)序錯(cuò)誤率設(shè)計(jì)映射所用的參數(shù)，具體的方法和步驟11相同；特別地，可以將測(cè)序序列與參考基因組之間的錯(cuò)誤匹配個(gè)數(shù)上界m調(diào)??；

步驟32、按照步驟31設(shè)計(jì)的參數(shù)，將測(cè)序序列映射至組裝基因組，將所有映射到多個(gè)位置的測(cè)序序列去除，得到單映測(cè)序序列形成的疊陣集；

步驟33、通過(guò)下述兩種方法之一獲得當(dāng)前預(yù)組裝疊陣集：

方法一：按照和步驟15相同的方法切割相對(duì)于當(dāng)前組裝基因組的單映測(cè)序序列疊陣；

方法二：找到組裝基因組上相鄰兩個(gè)連續(xù)為n的片段，截取位于它們中間的不包含n的堿基片段，以該堿基片段和單映射至該堿基片段的測(cè)序序列作為需要被延拓的疊陣，這些疊陣構(gòu)成當(dāng)前預(yù)組裝疊陣集；

步驟34、基于當(dāng)前預(yù)組裝疊陣集，利用單映射測(cè)序序列在當(dāng)前預(yù)組裝疊陣集上的坐標(biāo)信息，執(zhí)行步驟2實(shí)現(xiàn)迭代組裝。

在經(jīng)過(guò)一定次數(shù)的迭代組裝之后，輸出所得組裝基因組，作為方法的最終輸出。

以上為本發(fā)明的基本步驟。如按下述方案進(jìn)行修改，本發(fā)明在組裝基因組的同時(shí)，還可以構(gòu)建出組裝基因組上高雜合區(qū)域的雙倍體序列，并且同時(shí)輸出組裝基因組、雙倍體序列、雙倍體序列與組裝基因組的位置關(guān)系信息。所述修改方案為：以下述步驟a代替步驟226；以下述步驟b代替步驟23；在步驟3之后增加步驟c。

步驟a、整合x和y，將x向右延拓，推斷一致序列。完成此步驟的操作后可以得到一條一致序列，也可以得到兩條一致序列。如果得到兩條一致序列，每一條對(duì)應(yīng)于一個(gè)倍型。圖9示出了步驟a的步驟以及主要符號(hào)的含義。步驟a具體操作包括以下步驟：

如果y中含有來(lái)自x的測(cè)序序列，則進(jìn)入步驟a1，否則進(jìn)入步驟a6。

步驟a1、尋找一個(gè)同時(shí)屬于x和y的測(cè)序序列，記它在x和y中的坐標(biāo)為coorx和coory；對(duì)于y中的每一個(gè)測(cè)序序列，將其在y中的坐標(biāo)與(coorx-coory)相加，作為其在延拓后的疊陣x中的坐標(biāo)；對(duì)于同時(shí)包含于集合set以及x中的測(cè)序序列，如果其不在y中，則將其從x中去掉；

步驟a2、推斷x的一致序列，記作c(x)；

步驟a3、將y中的測(cè)序序列從g中去除，如果g還有剩余測(cè)序序列，且其中至少有一條與x中的某個(gè)測(cè)序序列存在有向邊，則進(jìn)入步驟a4；否則認(rèn)為x經(jīng)延拓后得到一條一致序列，結(jié)束步驟a的操作；

步驟a4、利用g中剩余的測(cè)序序列構(gòu)造一個(gè)疊陣hxr；記這些測(cè)序序列為r1，...，rm，一個(gè)優(yōu)選的構(gòu)造方式如下：

步驟a41、設(shè)定整數(shù)k1，...，km，將其初始化為0；將hxr設(shè)定為空集；

步驟a42、遍歷測(cè)序序列r1，...，rm，對(duì)于測(cè)序序列ri(1≤i≤m)，若ri與x中的一條測(cè)序序列v之間存在一條有向邊，則將ki更新為optimal_shift(ri，v)與v在x中的坐標(biāo)之和，并將ri加入hxr中，將其在hxr中的坐標(biāo)初始化為1；否則若ri與已存在于hxr中的一條測(cè)序序列rj之間存在一條有向邊，則將ki更新為kj+optimal_shift(ri，rj)，并將ri加入hxr中，將其在hxr中的坐標(biāo)初始化為1；

步驟a43、重復(fù)執(zhí)行步驟a42直到hxr中不再加入新的測(cè)序序列；

步驟a44、將k1，...，km中未被更新的值去掉，在余下的整數(shù)中找到最小值，記作kmin；對(duì)于測(cè)序序列ri(1≤i≤m)，如果ri位于kxr中，則將ri在hxr中的坐標(biāo)更新為ki-kmin+1；

步驟a5、推斷hxr的一致序列，記作c(hxr)；記錄數(shù)值kmin，作為hxr左端相對(duì)于x的平移量，也作為c(x)和c(hxr)的位置關(guān)系信息；將c(x)的第kmin個(gè)堿基標(biāo)記為分叉點(diǎn)，記之以符號(hào)pxr；結(jié)束步驟a的操作。

如果y中無(wú)來(lái)自x的測(cè)序序列，則執(zhí)行步驟a6：

步驟a6、分別推斷x和y的一致序列；如果y最左端的測(cè)序序列相對(duì)于x的先驗(yàn)坐標(biāo)值大于x的長(zhǎng)度，則將y的一致序列接在x的一致序列右邊，中間以若干個(gè)字符n相隔，形成一條一致序列，將其記為c(x)；否則保留兩條一致序列，將x的一致序列記為c(x)，將y的一致序列記為c(hxr)，同時(shí)記錄位于y最左端的測(cè)序序列的先驗(yàn)坐標(biāo)值priory_left，作為兩條一致序列的位置關(guān)系信息，將c(x)的第priory_left個(gè)堿基標(biāo)記為分叉點(diǎn)，記之以符號(hào)pxr。

步驟b、按照步驟24中描述的連接兩條序列的方案對(duì)相鄰的每?jī)蓚€(gè)疊陣的一致序列進(jìn)行連接，得到組裝基因組、連接區(qū)域的雙倍體序列，以及二者的位置關(guān)系信息，并將其輸出。

以下結(jié)合圖10(a)至圖10(d)，以兩個(gè)相鄰的疊陣a和b為例，闡述一致序列的連接，以及三項(xiàng)結(jié)果的輸出方法。以下分情況進(jìn)行闡述：

圖10(a)所示，a右端和b左端均存在雙倍型：判斷[c(a)]r和[c(b)]l、[c(a)]r和c(hbl)、c(har)和[c(b)]l、c(har)和c(hbl)能否成功連接；如果存在可以成功連接的序列對(duì)，則將連接部分長(zhǎng)度最大的序列對(duì)進(jìn)行連接，并將余下的兩個(gè)序列進(jìn)行連接；如果不存在可以成功連接的序列對(duì)，則任取一個(gè)序列對(duì)用n連接，將余下的兩個(gè)序列用n連接；如圖10(a)所示，被連接得到的兩條序列分別為[c(a)]rο[c(b)]l和c(har)οc(hbl)；；連續(xù)輸出[c(a)]c、[c(a)]rο[c(b)]l和[c(b)]c的每一個(gè)堿基作為組裝基因組；輸出c(har)οc(hbl)作為雙倍體序列；計(jì)算par和pbl在組裝基因組上的坐標(biāo)并輸出，作為雙倍體序列與組裝基因組的位置關(guān)系信息。

圖10(b)所示，a右端存在雙倍型，b左端不存在雙倍型：連接[c(a)]r和[c(b)]c得到[c(a)]rο[c(b)]c(或者連接c(har)和[c(b)]c得到c(har)ο[c(b)]c)，將[c(a)]rο[c(b)]c(或者c(har)ο[c(b)]c)輸出作為組裝基因組；輸出c(har)(或者[c(a)]r)作為雙倍體序列；計(jì)算par在組裝基因組上的坐標(biāo)并輸出，作為雙倍體序列與組裝基因組的位置關(guān)系信息。

圖10(c)所示，a右端不存在雙倍型，b左端存在雙倍型：連接[c(a)]c和[c(b)]l得到[c(a)]cο[c(b)]l(或者連接[c(a)]c和c(hbl)得到[c(a)]cοc(hbl))，將[c(a)]cο[c(b)]l(或者[c(a)]cοc(hbl))輸出作為組裝基因組；輸出c(hbl)(或者[c(b)]l)作為雙倍體序列；計(jì)算pbl在組裝基因組上的坐標(biāo)并輸出，作為雙倍體序列與組裝基因組的位置關(guān)系信息。

圖10(d)所示，a右端和b左端均不存在雙倍型：連接[c(a)]c和[c(b)]c得到[c(a)]cο[c(b)]c，輸出[c(a)]cο[c(b)]c作為組裝基因組。

從左至右順次掃描組裝基因組架構(gòu)，對(duì)每?jī)蓚€(gè)相鄰、且已被延拓的疊陣的一致序列進(jìn)行上述操作，即可得到組裝基因組，雙倍體序列，以及二者的位置關(guān)系信息。

步驟c、利用映射和基于序列兩兩比對(duì)的組裝方法構(gòu)建非連接處的雙倍體序列，作為步驟b所構(gòu)建的雙倍體序列的重要補(bǔ)充。圖11為步驟c的操作示意圖，主要操作包括測(cè)序序列再映射、局部化組裝、推斷局部化組裝所得疊陣的一致序列得到雙倍體序列。步驟c具體包括以下步驟：

步驟c1、將全體測(cè)序序列映射至步驟2得到的組裝基因組，和所有已構(gòu)建出的雙倍體序列。在映射所用的參數(shù)中，錯(cuò)誤匹配數(shù)目上界不應(yīng)設(shè)得過(guò)大，可以設(shè)為測(cè)序序列長(zhǎng)度的5％至6％。

步驟c2、對(duì)于步驟c1映射失敗的序列，再映射至組裝基因組；進(jìn)行再映射時(shí)，將錯(cuò)誤匹配數(shù)目的上界適度調(diào)大。如果進(jìn)行再映射后，映射失敗的序列比例仍然很大，可以將錯(cuò)誤匹配數(shù)目的上界繼續(xù)調(diào)大，對(duì)映射失敗的序列重復(fù)再映射。

步驟c3、對(duì)于任意一條在步驟c2中映射成功的測(cè)序序列r，進(jìn)行以下操作，以決定保留或者舍棄其映射結(jié)果：

步驟c31、如果r的同伴序列r′在步驟c1或者步驟c2中被成功映射至組裝基因組(而非步驟b構(gòu)建的雙倍體序列)，則進(jìn)入步驟c31；否則舍棄r的映射結(jié)果，結(jié)束步驟c3的操作。

步驟c32、記組裝基因組上r被成功映射到的位置為posi(1≤i≤s)，s表示r被成功映射到的位置數(shù)目；記組裝基因組上r′被成功映射到的位置數(shù)目為pos′j(1≤j≤s′)，s′表示r′被成功映射到的位置數(shù)目。如果存在唯一的地址對(duì)滿足：

(1)在和這兩個(gè)位置上，r與r′方向相反；

(2)和二者距離與r所屬片段庫(kù)的庫(kù)長(zhǎng)相近；

則保留r在的映射結(jié)果，舍棄余下s-1個(gè)映射結(jié)果；否則將r的全部s個(gè)映射結(jié)果舍棄。

步驟c4、對(duì)經(jīng)過(guò)步驟c3的篩選后被保留的測(cè)序序列按照在組裝基因組上的坐標(biāo)進(jìn)行分類。一個(gè)可行的分類方法如下：將所有被保留的測(cè)序序列按照其在組裝基因組上的坐標(biāo)從小到大進(jìn)行排序；將坐標(biāo)最小的測(cè)序序列分入第一類；按坐標(biāo)從小到大的順序依次掃描排好序的所有測(cè)序序列，如果一個(gè)測(cè)序序列的坐標(biāo)與上一個(gè)被掃描的測(cè)序序列坐標(biāo)之差小于設(shè)定的閾值，則將其與上一個(gè)被掃描的測(cè)序序列分入同一類，否則將其分入新的類中。

步驟c5、對(duì)于步驟c4得到的每一類，將其中所有測(cè)序序列按步驟232、步驟233、步驟234、步驟235所述方法構(gòu)建一個(gè)疊陣。構(gòu)建疊陣時(shí)所用到的先驗(yàn)序關(guān)系可以根據(jù)測(cè)序序列在組裝基因組上的坐標(biāo)獲得。推斷所構(gòu)建的疊陣的一致序列，作為一條雙倍體序列；記錄位于所構(gòu)建疊陣最左端與最右端的測(cè)序序列在組裝基因組上的坐標(biāo)，作為該雙倍體序列與組裝基因組的位置關(guān)系信息。

執(zhí)行完步驟c后，輸出以下內(nèi)容，作為方法的最終輸出包括：組裝基因組；步驟b和步驟c構(gòu)建的所有雙倍體序列；步驟b和步驟c構(gòu)建的每一條雙倍體序列與組裝基因組的位置關(guān)系信息。

本發(fā)明提出的方法適用基因組組裝。就基因組而言，本方法可以用于高雜合度、高重復(fù)度的基因組組裝；就測(cè)序序列而言，本方法適用于高通量、雙末端測(cè)序序列的組裝。如果存在與被測(cè)基因組相近的參考基因組，本方法可以通過(guò)映射、局部化組裝等步驟直接實(shí)現(xiàn)基因組組裝，以及雙倍體序列的構(gòu)建；在操作時(shí)，可以選擇一個(gè)參考基因組進(jìn)行上述操作，也可以選擇多個(gè)參考基因組，利用每一個(gè)參考基因組進(jìn)行上述操作，最后將結(jié)果整合，得到組裝基因組。如果不能找到合適的參考基因組，而且使用其他組裝方法得到的組裝基因組效果不理想時(shí)，可以將其他方法給出的組裝基因組作為參考基因組，利用本方法對(duì)被測(cè)基因組進(jìn)行重新組裝，進(jìn)行矯正。

本發(fā)明提出的方法還適用于不同基因組之間的比較。通過(guò)對(duì)疊陣進(jìn)行重新排列以及后續(xù)的局部化組裝，可以獲得被測(cè)基因組與參考基因組在結(jié)構(gòu)層面的差異信息，找到大片段的結(jié)構(gòu)變異，以及兩個(gè)物種基因組相似、或者變異密集的區(qū)段。

目前人類基因組已經(jīng)有若干版本的參考基因組，作為一個(gè)特別的應(yīng)用，本發(fā)明提出的方法可以基于人類參考基因組實(shí)現(xiàn)一個(gè)人的基因組組裝，并探測(cè)其與參考基因組之間存在的結(jié)構(gòu)層面的差異。這在個(gè)體化醫(yī)療中，對(duì)于實(shí)現(xiàn)基因組層面的疾病預(yù)測(cè)、診斷，可以起到重要的作用。

本發(fā)明提出的上述方法，具有以下優(yōu)點(diǎn)：

1、在設(shè)計(jì)思路上，本方法并不等同地看待所有測(cè)序序列，一次性完成組裝，而是先組裝出唯一性強(qiáng)、可靠性高的疊陣，通過(guò)這些疊陣獲得組裝基因組架構(gòu)，將重復(fù)區(qū)域和未組裝的部分局部化，之后進(jìn)行局部化的組裝，并構(gòu)建雙倍體序列；這樣有助于降低重復(fù)區(qū)域或者高雜合度給組裝帶來(lái)的不確定性；

2、通過(guò)將測(cè)序序列映射至參考基因組，根據(jù)映射結(jié)果得到一部分測(cè)序序列的疊落關(guān)系，用于實(shí)現(xiàn)被測(cè)基因組上與參考基因組相似度較高的地方的組裝，對(duì)組裝操作實(shí)現(xiàn)了部分程度的簡(jiǎn)化；

3、映射的參數(shù)設(shè)計(jì)有對(duì)應(yīng)的定量評(píng)估方法來(lái)指導(dǎo)，可以根據(jù)對(duì)被測(cè)基因組與參考基因組的差異預(yù)期、測(cè)序序列的長(zhǎng)度和質(zhì)量特征，定量地設(shè)計(jì)不同的參數(shù)，以實(shí)現(xiàn)靈敏度與特異度的平衡；

4、在映射后，僅保留映射到參考基因組上一個(gè)位置的測(cè)序序列，減少測(cè)序序列疊落關(guān)系中的不確定性；

5、測(cè)序序列映射、映射結(jié)果的過(guò)濾，對(duì)于不同測(cè)序序列都是獨(dú)立的，可以并行地實(shí)現(xiàn)；

6、在估計(jì)疊陣之間的距離時(shí)，同時(shí)利用不同片段庫(kù)的測(cè)序序列，增加估計(jì)時(shí)所用的樣本量；

7、利用參考基因組上具有自映射唯一性的區(qū)域過(guò)濾映射成功的測(cè)序序列，進(jìn)一步提高特異度，減少疊陣間的錯(cuò)誤連接；

8、排列疊陣后，被測(cè)基因組上的重復(fù)區(qū)域以及未被組裝的區(qū)域被局部化，表現(xiàn)為相鄰疊陣之間的空缺部分；已存在的疊陣可靠性較高，利用這些疊陣中的測(cè)序序列和雙末端信息，尋找屬于空缺部分的測(cè)序序列，使測(cè)序序列和空缺部分能夠更加準(zhǔn)確地對(duì)應(yīng)；

9、經(jīng)過(guò)局部化后，在延拓每一個(gè)疊陣時(shí)，所需組裝的測(cè)序序列數(shù)目大幅度減少，使得進(jìn)行延拓可以采用基于測(cè)序序列疊落關(guān)系的組裝方法進(jìn)行，從而保留了測(cè)序序列整體信息，無(wú)需將其切為k-mer；同時(shí)可以并行地對(duì)每一個(gè)疊陣進(jìn)行延拓，增加方案執(zhí)行的效率；

10、延拓疊陣、連接一致序列后，可以得到雙倍體序列，以及雙倍體和組裝基因組之間的位置關(guān)系信息；

11、構(gòu)建非連接處的雙倍體序列時(shí)，基于再映射和雙末端信息尋找屬于另一個(gè)倍型的測(cè)序序列；在映射時(shí)允許的錯(cuò)誤匹配數(shù)目較多，雙末端信息的應(yīng)用有助于糾正由此導(dǎo)致的錯(cuò)誤映射；另外，找到的屬于另一個(gè)倍型的測(cè)序序列具有在組裝基因組上的再映射坐標(biāo)，雙倍體序列與組裝基因組的位置關(guān)系可以較為準(zhǔn)確地獲得；

12、延拓疊陣，連接一致序列可以并行地實(shí)現(xiàn)；尋找非連接處的雙倍體序列，也可以將測(cè)序序列根據(jù)再映射坐標(biāo)進(jìn)行分類，并行地實(shí)現(xiàn)。

根據(jù)本發(fā)明的實(shí)施例，還提供了一種應(yīng)用上述組裝方法中任意一種組裝方法的變異探測(cè)方法。

根據(jù)本發(fā)明實(shí)施例的變異探測(cè)方法包括：

根據(jù)組裝方法對(duì)不同樣本的基因組之間的結(jié)構(gòu)變異情況進(jìn)行探測(cè)，探測(cè)的信息包括對(duì)單映射測(cè)序序列疊陣集進(jìn)行切割時(shí)形成的斷點(diǎn)信息。

根據(jù)本發(fā)明的實(shí)施例，還提供了一種基因組序列的組裝系統(tǒng)。

如圖12所示，根據(jù)本發(fā)明實(shí)施例的組裝系統(tǒng)包括：

映射模塊121，用于通過(guò)預(yù)定的映射算法將樣品的被測(cè)基因組的測(cè)序序列映射到參考基因組，得到單映射測(cè)序序列疊陣集，其中，樣品的測(cè)序序列為利用高通量測(cè)序技術(shù)測(cè)得，參考基因組已知并與樣品的基因組相近；

篩選模塊122，用于基于經(jīng)過(guò)預(yù)處理的參考基因組對(duì)單映射測(cè)序序列疊陣集中的測(cè)序序列進(jìn)行篩選，所得篩選結(jié)果根據(jù)覆蓋度再次篩選，得到篩選后的單映射測(cè)序序列疊陣集；

切割模塊123，用于通過(guò)單方向測(cè)序序列信息對(duì)篩選后的單映射測(cè)序序列疊陣集進(jìn)行切割，得到初始預(yù)組裝疊陣集；

架構(gòu)模塊124，用于確定初始預(yù)組裝疊陣集中每個(gè)疊陣的相對(duì)位置，形成組裝基因組架構(gòu)；

延拓模塊125，用于對(duì)組裝基因組架構(gòu)中的每個(gè)疊陣進(jìn)行延拓，得到每個(gè)疊陣的一致序列；

連接模塊126，用于將組裝基因組架構(gòu)中的相鄰疊陣的一致序列中符合預(yù)定連接規(guī)則的一致序列進(jìn)行連接，得到樣品的當(dāng)前的組裝基因組；

調(diào)整映射模塊127，用于根據(jù)被測(cè)基因組上同源序列的差異調(diào)整預(yù)定的映射算法的預(yù)定的映射參數(shù)，通過(guò)調(diào)整的該預(yù)定的映射算法將樣品的被測(cè)基因組的測(cè)序序列映射到當(dāng)前的組裝基因組，得到當(dāng)前預(yù)組裝疊陣集；

架構(gòu)模塊124，延拓模塊125和連接模塊126進(jìn)一步用于對(duì)調(diào)整映射模塊127中的當(dāng)前預(yù)組裝疊陣集進(jìn)行操作。

綜上所述，借助于本發(fā)明的上述技術(shù)方案，通過(guò)將被測(cè)基因組的測(cè)序序列和參考基因組進(jìn)行映射，并對(duì)映射結(jié)果進(jìn)行切割，以及將切割后的疊陣進(jìn)行組裝和延拓，從而實(shí)現(xiàn)了測(cè)序序列的高效連接，實(shí)現(xiàn)測(cè)序序列的基因組裝。

總之，本發(fā)明公開了一種基因組序列的組裝方法、和相應(yīng)的結(jié)構(gòu)變異探測(cè)方法，以及基因組組裝系統(tǒng)，該組裝方法包括：通過(guò)設(shè)計(jì)序列映射的唯一性準(zhǔn)則，將被測(cè)基因組的測(cè)序序列和參考基因組進(jìn)行映射，并對(duì)映射結(jié)果進(jìn)行恰當(dāng)?shù)那懈?，形成預(yù)組裝疊陣集。然后根據(jù)單映射序列在組裝疊陣集上的坐標(biāo)和同伴關(guān)系估計(jì)基因組的架構(gòu)，并根據(jù)組裝疊陣集上單映序列的坐標(biāo)和它們的同伴序列將疊陣向外延拓。延拓采用以下三部曲算法：1.序列兩兩比對(duì)；2.用圖論方法整合兩兩比對(duì)結(jié)果形成延拓后的疊陣；3.基于所述延拓后的疊陣定義延拓的一致序列。上述延拓對(duì)各個(gè)疊陣以并行方式計(jì)算執(zhí)行。延拓后的相鄰疊陣一致序列經(jīng)過(guò)比對(duì)判別，如果存在重疊，就將它們連接，從而完成一輪的基因組拼接。所得到的當(dāng)前組裝基因組作為下一輪的參考基因組，通過(guò)調(diào)整序列映射的唯一性準(zhǔn)則，重復(fù)以上拼接步驟，改進(jìn)基因組的組裝結(jié)果。所測(cè)基因組相對(duì)于參考基因組的結(jié)構(gòu)變異，在拼接的過(guò)程中同時(shí)被探測(cè)出來(lái)。

這項(xiàng)知識(shí)產(chǎn)權(quán)的研發(fā)得到了中國(guó)科學(xué)院b類先導(dǎo)專項(xiàng)“動(dòng)物復(fù)雜性狀的進(jìn)化解析與調(diào)控"課題xdb13040600的資助，得到了國(guó)家自然科學(xué)基金委員會(huì)重大研究計(jì)劃培育項(xiàng)目91530105、91130008的資助，以及中國(guó)科學(xué)院國(guó)家數(shù)學(xué)與交叉科學(xué)中心的各種支持。

以上所述僅為本發(fā)明的較佳實(shí)施例而已，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內(nèi)，所作的任何修改、等同替換、改進(jìn)等，均應(yīng)包含在本發(fā)明的保護(hù)范圍之內(nèi)。