雜合基因組處理方法
【專利摘要】本發(fā)明公開了一種雜合基因組處理方法，其包括采用全基因組鳥槍法對目標(biāo)物種進(jìn)行至少200倍的超高深度測序，以獲得有效的讀長短序列Reads；對所述有效的Reads執(zhí)行組裝和構(gòu)建支架序列Scaffold，以獲得帶有冗余序列的基因組圖譜；對所述帶有冗余序列的基因組圖譜執(zhí)行雜合識別處理，從而去除雜合區(qū)域中冗余的Scaffold并保留雜合區(qū)域信息以獲得全基因組圖譜。通過本發(fā)明的方法實(shí)現(xiàn)了快速組裝成完整高質(zhì)量的雜合基因組圖譜，縮短了繪制周期和成本。
【專利說明】雜合基因組處理方法

【技術(shù)領(lǐng)域】
[0001]本發(fā)明涉及基因工程【技術(shù)領(lǐng)域】和基因組生物信息學(xué)分析領(lǐng)域。特別地，本發(fā)明涉及雜合基因組處理，以獲得物種的全基因組序列圖譜。

【背景技術(shù)】
[0002]全基因組de novo測序也叫從頭測序，不需要任何參考序列即可對某個(gè)物種進(jìn)行測序，用生物信息學(xué)分析方法進(jìn)行拼接、組裝，從而獲得該物種的全基因組序列圖譜。獲得一個(gè)物種的全基因組序列是加快對此物種了解的重要捷徑。
[0003]隨著人類基因組計(jì)劃啟動(dòng)，世界邁出了基因組測序的重要一步；各種模式的動(dòng)植物線蟲、擬南芥、小鼠等紛紛進(jìn)行了基因組測序；由于技術(shù)水平限制，基因組測序還處于初級階段，在一些重要物種測序完成后，第一代基因組測序也慢慢走向成熟，伴隨而來的是利用第一代測序技術(shù)對各種重要物種基因組的測序熱潮，各種重要?jiǎng)又参锘蚪M相繼被測序。不論這些數(shù)字如何令人印象深刻，第一代測序技術(shù)在速度和成本方面都已達(dá)到了極限，因此，需要開發(fā)全新的技術(shù)來突破這些局限。另外，第一代測序技術(shù)也遇到了一個(gè)很大的瓶頸，那就是基因組雜合的問題?；蚪M雜合的情況廣泛存在于具有重要經(jīng)濟(jì)或科研價(jià)值的物種中，而傳統(tǒng)的雜合基因組圖譜繪制方案成本高、周期長，嚴(yán)重制約了這些物種基因組的解讀。


【發(fā)明內(nèi)容】

[0004]本發(fā)明要解決的技術(shù)問題在于，針對現(xiàn)有技術(shù)中傳統(tǒng)的雜合基因組圖譜繪制方案成本高、周期長，嚴(yán)重制約了物種基因組的解讀的缺陷，提供一種雜合基因組處理方法。
[0005]本發(fā)明解決其技術(shù)問題所采用的技術(shù)方案是:構(gòu)造一種雜合基因組處理方法，其包括:
采用全基因組鳥槍法對目標(biāo)物種進(jìn)行至少200倍的超高深度測序，以獲得有效的讀長短序列Reads ；
對所述有效的Reads執(zhí)行組裝和構(gòu)建支架序列Scaffold，以獲得帶有冗余序列的基因組圖譜；
對所述帶有冗余序列的基因組圖譜執(zhí)行雜合識別處理，從而去除雜合區(qū)域中冗余的Scaffold并保留雜合區(qū)域信息以獲得全基因組圖譜。
[0006]在本發(fā)明的雜合基因組處理方法中，所述雜合識別處理包括基于k-mer分布圖來識別所述雜合區(qū)域。
[0007]在本發(fā)明的雜合基因組處理方法中，所述基于k-mer分布圖來識別所述雜合區(qū)域包括:
通過k-mer分布圖，識別所述帶有冗余序列的基因組圖譜中的Unique k-mer區(qū)域,其中，所述Unique k_mer表示在基因組中只出現(xiàn)一次的k_mer,k_mer表示定長為指定值k的短序列，k為正整數(shù)；根據(jù)所述Unique k_mer區(qū)域來確定所述雜合區(qū)域。
[0008]在本發(fā)明的雜合基因組處理方法中，所述根據(jù)所述Unique k-mer區(qū)域來確定所述雜合區(qū)域包括:
統(tǒng)計(jì)每條Scaffold上的Unique k-mer的數(shù)量N以及那些在已統(tǒng)計(jì)過的一條Scaffold上出現(xiàn)過的Unique k-mer的數(shù)量M，其中N和M均為正整數(shù)；
如果M/N的值達(dá)到預(yù)定值，則將這兩條Scaffold確定為對應(yīng)的雜合區(qū)域。
[0009]在本發(fā)明的雜合基因組處理方法中，至少基于組裝的長度來去除雜合區(qū)域中冗余的 Scaffoldo
[0010]在本發(fā)明的雜合基因組處理方法中，在構(gòu)建支架序列Scaffold之前，還包括構(gòu)建重疊群序列Contig。
[0011]在本發(fā)明的雜合基因組處理方法中，還包括對所述全基因組圖譜進(jìn)行基因組注釋和/或基因組進(jìn)化分析。
[0012]在本發(fā)明的雜合基因組處理方法中，所述基因組注釋包括R印eat注釋、基因預(yù)測、基因功能注釋和ncRNA注釋；所述基因組進(jìn)化分析包括基因家族鑒定、物種系統(tǒng)發(fā)育樹構(gòu)建、物種分歧時(shí)間估算、基因組共線性分析和全基因組復(fù)制分析。
[0013]在本發(fā)明的雜合基因組處理方法中，采用全基因組鳥槍法對目標(biāo)物種進(jìn)行至少200倍的超高深度測序，以獲得有效的讀長短序列Reads包括:
構(gòu)建文庫；
采用全基因組鳥槍法對所構(gòu)建的文庫進(jìn)行至少200倍的超高深度測序，以獲得原始的Reads ；
對所述原始的Reads進(jìn)行過濾，以獲得有效的Reads。
[0014]在本發(fā)明的雜合基因組處理方法中，構(gòu)建文庫包括構(gòu)建小片段文庫和構(gòu)建大片段文庫，其中，小片段文庫包括長度不同的多個(gè)插入片段文庫，大片段文庫包括長度不同的多個(gè)配對Mate Pair文庫。
[0015]在本發(fā)明的雜合基因組處理方法中，構(gòu)建小片段文庫包括:
將檢測合格的目標(biāo)物種的DNA樣品打斷，以獲得DNA片段；
對所述DNA片段進(jìn)行末端修復(fù)；
將末端修復(fù)的DNA片段的3’端連接堿基A ；
對連接堿基A的DNA片段連接上接頭Adapter ；
根據(jù)目的片段的長度，對連接上接頭Adapter的DNA片段執(zhí)行PCR反應(yīng)和PCR擴(kuò)增，以構(gòu)建小片段文庫。
[0016]在本發(fā)明的雜合基因組處理方法中，構(gòu)建大片段文庫包括:
將目標(biāo)物種的基因組DNA打斷，以獲得DNA片段； 對所述DNA片段進(jìn)行末端修復(fù)和生物素標(biāo)記；
對經(jīng)末端修復(fù)和生物素標(biāo)記的DNA片段進(jìn)行片段選擇，以獲得目的DNA片段；
對目的DNA片段進(jìn)行環(huán)化并消化線性DNA ；
將環(huán)化的目的DNA片段打斷成預(yù)定長度的片段，并捕獲帶有生物素標(biāo)記的片段； 對捕獲的片段進(jìn)行末端修飾和連接上特定的接頭Adapter，以構(gòu)建大片段文庫。
[0017]實(shí)施本發(fā)明的技術(shù)方案，本發(fā)明不需要任何參考序列下，基于新一代高通量測序技術(shù),通過對目標(biāo)物種進(jìn)行超高深度測序(至少200倍)；不需要構(gòu)建BAC文庫，成本大大降低；本發(fā)明的技術(shù)方案只需要6個(gè)月，而現(xiàn)有的技術(shù)方案需要I年，因此周期大大縮短；將極大推動(dòng)雜合基因組研究的進(jìn)展，為快速獲得物種的高質(zhì)量基因組圖譜提供了強(qiáng)有力保障。

【專利附圖】

【附圖說明】
[0018]下面將結(jié)合附圖及實(shí)施例對本發(fā)明作進(jìn)一步說明，附圖中:
圖1是根據(jù)本發(fā)明的實(shí)施例的雜合基因組處理方法的示意流程圖；
圖2是根據(jù)本發(fā)明的實(shí)施例的雜合識別處理的示意圖；
圖3是根據(jù)本發(fā)明的另一實(shí)施例的雜合基因組處理方法的示意圖；
圖4是根據(jù)本發(fā)明的實(shí)施例的構(gòu)建小片段文庫的流程圖；
圖5是根據(jù)本發(fā)明的實(shí)施例的構(gòu)建大片段文庫的流程圖；
圖6是根據(jù)本發(fā)明的實(shí)施例的雜合基因組處理方法中信息分析流程的示意圖；
圖7是根據(jù)本發(fā)明的實(shí)施例的17-mer深度頻率分布圖；
圖8是根據(jù)本發(fā)明的實(shí)施例的GC含量與測序深度關(guān)系分布圖；
圖9是根據(jù)本發(fā)明的實(shí)施例的四種物種的基因組GC含量分布的示意圖；
圖10是根據(jù)本發(fā)明的實(shí)施例的測序深度分布圖；
圖11是根據(jù)本發(fā)明的實(shí)施例的多個(gè)物種的蛋白家族的同源性的示意圖；
圖12是根據(jù)本發(fā)明的實(shí)施例的系統(tǒng)發(fā)育樹的示意圖；
圖13是根據(jù)本發(fā)明的實(shí)施例的分化時(shí)間和替換速率的估算的示意圖；
圖14是根據(jù)本發(fā)明的實(shí)施例的大熊貓、狗和人三個(gè)基因組的保守序列長度關(guān)系的示意圖；
圖15是根據(jù)本發(fā)明的實(shí)施例的兩個(gè)物種的基因組共線性的示意圖；
圖16是根據(jù)本發(fā)明的實(shí)施例的兩個(gè)物種的共線性示意圖；
圖17是根據(jù)本發(fā)明的實(shí)施例的大片段復(fù)制長度、相似性的分布圖；
圖18是根據(jù)本發(fā)明的實(shí)施例的基因組自身共線性比對分析結(jié)果的示意圖；
圖19是根據(jù)本發(fā)明的實(shí)施例的基因組自身共線性比對分析結(jié)果的另一示意圖；
圖20是根據(jù)本發(fā)明的實(shí)施例的兩個(gè)物種的4DTv分布圖。

【具體實(shí)施方式】
[0019]本發(fā)明的目的、特征和特性以及結(jié)構(gòu)的有關(guān)元件的操作方法和功能和部件的組合在參照附圖閱讀下列描述和以上權(quán)利要求時(shí)將變得更明顯。要明確地理解附圖僅僅是為了說明和描述的目的并且不意在作為本發(fā)明的限制。
[0020]在本發(fā)明的雜合基因組處理方法中，不需要任何參考序列下，基于新一代高通量測序技術(shù)，通過對目標(biāo)物種進(jìn)行超高深度測序(一般在200倍以上)，然后采用新版S0APdenovo2組裝軟件，結(jié)合雜合識別處理，從而獲得目標(biāo)物種的全基因組序列圖譜。本發(fā)明的實(shí)施例可以包括全基因組鳥槍法(WGS)測序、基因組組裝、基因組注釋和進(jìn)化分析，以及個(gè)性化分析，旨在解決對于雜合基因組的難題。
[0021]如圖1所示，在本發(fā)明的實(shí)施例中， 雜合基因組處理方法包括:5101:采用全基因組鳥槍法對目標(biāo)物種進(jìn)行至少200倍的超高深度測序，以獲得有效的讀長短序列Reads ；
5102:對所述有效的Reads執(zhí)行組裝和構(gòu)建支架序列Scaffold，以獲得帶有冗余序列的基因組圖譜；
5103:對所述帶有冗余序列的基因組圖譜執(zhí)行雜合識別處理，從而去除雜合區(qū)域中冗余的Scaffold并保留雜合區(qū)域信息以獲得全基因組圖譜。在特定實(shí)施例中，如圖2所示，雜合識別處理過程可以包括:通過k-mer (其表示定長為指定值k的短序列)分布圖確定Unique k_mer,即在基因組中只出現(xiàn)一次的k_mer,其中，Unique k_mer來自序列中的Unique區(qū)域,其代表著獨(dú)一無二的基因信息，對拼接有著重要意義，而Repeat k-mer可能來自序列中的重復(fù)區(qū)域，這些重復(fù)區(qū)域可能帶來潛在的錯(cuò)拼。接著，將WGS測序的有效的Reads打斷成k-mer,統(tǒng)計(jì)k_mer的頻率并繪制曲線(在一個(gè)實(shí)施例中,k_mer曲線繪制的時(shí)候可以采用17-mer),利用k-mer分布查找基因組中的Unique k-mer。從概率上來講，> 99%的Unique k-mer分布在1.5倍主峰之前的區(qū)域,所以把此區(qū)域內(nèi)的k_mer定義為Uniquek-mer (如圖2所示矩形框內(nèi)部為Unique K-mer);將Scaffold序列從大到小排列,然后將Unique k-mer比對回Scaffold序列，Unique k-mer富集的區(qū)域稱為Unique區(qū)域，例如圖2中的A、A’。而圖2中B或C僅為示例基因組存在重復(fù)序列等的其他類型序列(repeatk-mer較多的區(qū)域),并不是通過序列比對得到這些區(qū)域信息；統(tǒng)計(jì)每條Scaffold上Uniquek-mer的數(shù)量N,以及那些在已統(tǒng)計(jì)過的Scaffold上出現(xiàn)過的Unique k-mer的數(shù)量Μ,如果Μ/Ν的值達(dá)到預(yù)定值(在此，該預(yù)定值可根據(jù)實(shí)際實(shí)現(xiàn)和需要而靈活設(shè)置，為此不進(jìn)行特別限定)，說明這條序列上的Unique區(qū)域與已統(tǒng)計(jì)過的其他較長Scaffold的Unique區(qū)域存在雜合現(xiàn)象。如果兩條Scaffold被確定為對應(yīng)雜合區(qū)域，將至少基于組裝長度，將冗余的Scaffold去除，例如可以將一條組裝質(zhì)量較差，組裝長度的短的Scaffold視為冗余，從基因組圖譜中去除，保留在另一份文件中，以最終組裝成完整高質(zhì)量基因組圖譜，同時(shí)保留雜合區(qū)域信息，以便后續(xù)進(jìn)行基因組注釋和基因組進(jìn)化分析，以及個(gè)性化分析。在一個(gè)實(shí)施例中，基因組注釋包括R印eat注釋、基因預(yù)測、基因功能注釋和ncRNA注釋；基因組進(jìn)化分析包括基因家族鑒定、物種系統(tǒng)發(fā)育樹構(gòu)建、物種分歧時(shí)間估算、基因組共線性分析和全基因組復(fù)制分析。由此，本發(fā)明的技術(shù)方案與WGS+BAC to BAC方案相比不需要構(gòu)建BAC文庫，成本大大降低，以及只需要6個(gè)月，而WGS+BAC to BAC方案需要I年。
[0022]另外，在一實(shí)施例中，在構(gòu)建支架序列Scaffold之前，還包括構(gòu)建重疊群序列Contig。如圖3所示，根據(jù)本發(fā)明的另一實(shí)施例的雜合基因組處理方法包括獲取目標(biāo)物種的基因組DNA、針對該基因組DNA構(gòu)建WGS測序文庫、進(jìn)行至少200倍超高深度測序、構(gòu)建Contig、構(gòu)建Scaffold以獲得帶冗余序列的基因組圖譜、識別Unique k_mer、雜合識別處理以獲得全基因組圖譜。
[0023]在一個(gè)實(shí)施例中，采用全基因組鳥槍法對目標(biāo)物種進(jìn)行至少200倍的超高深度測序，以獲得有效的讀長短序列Reads包括:構(gòu)建文庫；采用全基因組鳥槍法對所構(gòu)建的文庫進(jìn)行至少200倍的超高深度測序，以獲得原始的Reads ;對所述原始的Reads進(jìn)行過濾，以獲得有效的Reads。
[0024]在一個(gè)實(shí)施例中，構(gòu)建文庫包括構(gòu)建小片段文庫和構(gòu)建大片段文庫，其中，小片段文庫包括長度不同的多個(gè)插入片段文庫，大片段文庫包括長度不同的多個(gè)Mate Pair文庫。在本發(fā)明的實(shí)施例中，小片段數(shù)據(jù)指插入片段長度小于I kb，大片段數(shù)據(jù)指插入片段長度大于I kb ο
[0025]在特定的實(shí)施例中，構(gòu)建17(^口、50(^口、80(^口、21^、51^、101^、201^和 40kb 共 8 種 WGS
測序文庫。根據(jù)基因組大小和重復(fù)序列的具體特點(diǎn)，制定詳細(xì)的基因組測序策略，整體測序深度在200倍以上，以保證單堿基的準(zhǔn)確性和基因組的完整性。小片段文庫可包括170 bp、500 bp >800 bp插入片段文庫；大片段文庫可包括2 kb>5 kb、10 kb >20 kb、40kb配對(MatePair)文庫。
[0026]如圖4所示，根據(jù)本發(fā)明的實(shí)施例的構(gòu)建小片段文庫包括:
S401:將檢測合格的目標(biāo)物種的DNA樣品打斷(即對DNA樣品進(jìn)行片段化)，以獲得DNA片段；其中，可以利用Nanodrop、Gel_Electrophotometric等對目標(biāo)物種的DNA樣品進(jìn)行檢測，經(jīng)檢測合格的DNA樣品用于文庫制備。另外，可以使用超聲打斷DNA樣品。
[0027]S402:對DNA片段進(jìn)行末端修復(fù)。在一個(gè)實(shí)施例中，可以使用T4 DNA聚合酶和Klenow聚合酶修復(fù)末端，制備平末端。
[0028]S403:將末端修復(fù)的DNA片段的3’端連接堿基A。在一個(gè)實(shí)施例中，在聚合酶體系的作用下，使上一步得到的修復(fù)產(chǎn)物在3 ’端加上堿基A,為下一步的接頭Adapter連接做準(zhǔn)備。
[0029]S404:對連接堿基A的DNA片段連接上接頭Adapte。在一個(gè)實(shí)施例中，配置T4 DNA連接酶反應(yīng)體系，在Thermomixer中適溫反應(yīng)一定時(shí)間使adapter和加“A”產(chǎn)物連接。
[0030]S405:DNA片段大小選擇。在一個(gè)實(shí)施例中，根據(jù)目的片段大小的不同制備相應(yīng)的瓊脂糖凝膠，對上一步的連接產(chǎn)物進(jìn)行PCR反應(yīng)。
[0031]S406 =PCR擴(kuò)增及產(chǎn)物純化。在一個(gè)實(shí)施例中，配置PCR反應(yīng)體系對上一步的切膠回收產(chǎn)物進(jìn)行擴(kuò)增，其中PCR酶選用Phus1n DNA Polymerase,瓊脂糖凝膠純化,從而完成小片段文庫構(gòu)建。
[0032]如圖5所示，根據(jù)本發(fā)明的構(gòu)建大片段文庫包括:
S501:將目標(biāo)物種的基因組DNA隨機(jī)打斷到特定大小，以得到DNA片段。例如，在一個(gè)實(shí)施例中，可選擇2-10 kb的范圍。
[0033]S502:對DNA片段進(jìn)行末端修復(fù)和生物素標(biāo)記。在一個(gè)實(shí)施例中，T4 DNA聚合酶可將3’突出末端平滑化，將5’突出末端補(bǔ)平。同時(shí)通過置換合成法標(biāo)記探針，T4 PNK可將DNA5’端磷酸化，以便進(jìn)行連接反應(yīng)，并進(jìn)行DNA末端標(biāo)記，用T4多聚核酸激酶除去3’磷酸基團(tuán)。Klenow DNA聚合酶可補(bǔ)平5’突出末端，形成平末端；切除3’突出末端，形成平末端。利用這三種酶的作用，可達(dá)到DNA片段修復(fù)和生物素標(biāo)記的目的。
[0034] S503 =DNA片段大小選擇，以獲得目的DNA片段。
[0035]S504:對目的DNA片段進(jìn)行環(huán)化。在一個(gè)實(shí)施例中，取適量回收后的DNA進(jìn)行環(huán)化的酶反應(yīng)；
S505:線形DNA的消化，其中，將環(huán)化后的DNA樣品加入消化反應(yīng)體系中，消化線性
DNA。
[0036]S506:將環(huán)化的目的DNA片段打斷成預(yù)定長度的片段，并捕獲帶有生物素標(biāo)記的片段。在一個(gè)實(shí)施例中，破碎環(huán)化的DNA并純化:把環(huán)化后的DNA分子打斷成400-600 bp的片段，并通過帶有鏈親和霉素的磁珠把那些帶有生物素標(biāo)記的片段捕獲。
[0037]S507:對捕獲的片段進(jìn)行末端修飾和連接上特定的接頭Adapter，以構(gòu)建大片段文庫。
[0038]在一個(gè)實(shí)施例中，對小片段文庫和大片段文庫的整體測序深度均至少為200倍。小片段文庫測序策略為101 PE ;大片段文庫測序策略為50 PE或91 PE0采取IlluminaHiseq2000測序，可得到很高的基因組單堿基測序深度和基因組物理覆蓋度，從而帶來高的單堿基準(zhǔn)確率，并確保基因組的完整性。
[0039]然而，測序得到的原始Reads，并不都是有效的，里面含有帶接頭的，重復(fù)的，測序質(zhì)量很低的Reads,這些Reads會(huì)影響組裝和后續(xù)分析,必須對原始的Reads過濾,得到有效Reads。過濾的Reads方法如下:
a)截去5’端和3’端堿基GC分叉嚴(yán)重的部分；
b)去除N堿基(含量超過2%)或者是polyA的Reads ；
c)對于小片段數(shù)據(jù)去除低質(zhì)量(質(zhì)量值〈=7)堿基數(shù)目達(dá)到一定程度的Reads(取不大于40%);對于大片段數(shù)據(jù)去除低質(zhì)量堿基數(shù)目達(dá)到一定程度的Reads (取不大于60%)；。
[0040]d)去除有Adapter污染的Reads (與Adapter序列至少10 bp比對上，且失配數(shù)不多于3個(gè)。
[0041]e)去除Readl和Read2有重疊的Reads (Readl和Read2重疊至少10 bp，且失配低于10%。對測通的Reads不進(jìn)行此操作。
[0042]f )去除重復(fù)的Reads (Readl和Read2完全一樣才算是重復(fù))。
[0043]利用k-mer進(jìn)行數(shù)據(jù)糾錯(cuò)，對于一個(gè)長度為L的序列可以得到(L_k+1)個(gè)長度為k的序列。在數(shù)據(jù)量足夠的情況下，從Reads中逐堿基取出的所有k-mer能夠遍歷整個(gè)基因組，并且k-mer出現(xiàn)的頻數(shù)是服從泊松分布的。測序錯(cuò)誤會(huì)導(dǎo)致新的k-mer出現(xiàn)，一般來說這些新的k-mer頻數(shù)都是比較低的。當(dāng)測序量足夠大的情況下，可以認(rèn)為頻數(shù)低的k-mer是因?yàn)闇y序錯(cuò)誤導(dǎo)致的。因?yàn)榇笃卧诮◣斓臅r(shí)候經(jīng)過的環(huán)化的過程，而且在組裝過程中只起到定位的作用，所以糾錯(cuò)僅針對小片段數(shù)據(jù)進(jìn)行。先利用過濾后數(shù)據(jù)建立k-mer頻數(shù)表，設(shè)置切割將k-mer分為高頻的和低頻的。對于有低頻k-mer出現(xiàn)的Reads，通過改變某些堿基可以使得整個(gè)Reads上的k-mer都為高頻，那么就認(rèn)為糾正了一些測序?qū)е碌腻e(cuò)誤。完成基因組組裝后，為了評價(jià)單堿基覆蓋深度并反映堿基準(zhǔn)確性，進(jìn)行堿基深度分析。采用比對程序SOAPaligner將過濾之后的Reads比對回拼接的基因組序列上，然后根據(jù)比對結(jié)果統(tǒng)計(jì)每個(gè)堿基被覆蓋的次數(shù)，從而可以得到各種測序深度的堿基占全基因組的百分比。華大基因自主研發(fā)了針對于Illumina Hiseq2000技術(shù)特點(diǎn)的短序列組裝分析軟件(SOAP軟件)，并已成功應(yīng)用于黃瓜、熊貓、白菜、螞蟻、馬鈴薯、谷子、牦牛、牡蠣等上百個(gè)基因組的組裝，積累了豐富的基因組組裝經(jīng)驗(yàn)，能高效準(zhǔn)確地對HiSeq2000測序得到的海量數(shù)據(jù)進(jìn)行拼接組裝。目前Soapdenovo已經(jīng)升級到Soapdenovo2,較之前的版本有多項(xiàng)功能改進(jìn)，組裝效果更好。在Contig構(gòu)建中，Soapdenovo2具有可選的多k-mer的功能,多k-mer是soapdenovo2的新特性。可以在Contig構(gòu)建中使用梯度k_mer,較之前版本的單一 k_mer組裝，效果更好。
[0044]如圖6所示，在根據(jù)本發(fā)明的實(shí)施例的雜合基因組處理方法中信息分析流程中涉及對大片段文庫和小片段文庫進(jìn)行測序所得的原始Reads進(jìn)行數(shù)據(jù)過濾、基因組組裝、基因組注釋和基因組進(jìn)化分析?；蚪M組裝涉及采用S0APdenOVO2組裝軟件進(jìn)行組裝、GC-深度分析、GC含量分析、測序深度分析、基因組組裝結(jié)果評估(常染色體區(qū)域覆蓋度評估和基因區(qū)覆蓋度評估)。基因組注釋涉及R印eat注釋、基因預(yù)測、基因功能注釋和ncRNA注釋。基因組進(jìn)化分析涉及基因家族鑒定、物種系統(tǒng)發(fā)育樹構(gòu)建、物種分歧時(shí)間估算、基因組共線性分析和全基因組復(fù)制分析。
[0045]在一個(gè)實(shí)施例中，如圖7所示，對于GC-深度分析，可執(zhí)行k-mer分析估計(jì)基因組大小。具體地，從一段連續(xù)序列中迭代地選取長度為k個(gè)堿基的序列，若Reads長度為L，k-mer長度為k,那么可以得到L_k+1個(gè)k-mer。例如,可取k為17來進(jìn)行分析。分析結(jié)果如圖7:圖中橫坐標(biāo)為深度，縱坐標(biāo)為各深度下的k-mer種類占所有k-mer種類的比例，即頻率。
[0046]使用高質(zhì)量測序數(shù)據(jù)，逐堿基取17-mer獲得深度頻數(shù)分布圖，如圖7示出深度峰值的位置，因此，根據(jù)公式Genome Size = k_mer_num/Peak _ckpth,可以估算出該物種的基因組大小。本發(fā)明的技術(shù)方案與現(xiàn)有技術(shù)方案最大的區(qū)別是k-mer大小的設(shè)置和一些參數(shù)的調(diào)節(jié):在組裝的時(shí)候用的是大k-mer進(jìn)行組裝，一般是取k=61以上，而現(xiàn)有技術(shù)大多取的是k=30多或40多。大k-mer能跨過更長的重復(fù)序列和更多的雜合位點(diǎn)，但是它對數(shù)據(jù)量的要求非常高，以61-mer為例，它至少需要150-160倍以上的測序深度，而63_mer則要求160-180倍測序深度。k值取的越大，要求的測序量越大。
[0047]在一個(gè)實(shí)施例中，如圖8所示，對于GC含量分析，各個(gè)物種基因組序列中GC堿基或AT堿基的分布呈現(xiàn)不同特征，而且不同的GC含量可能影響測序隨機(jī)性。為了分析該基因組序列的堿基分布特征和測序隨機(jī)性，進(jìn)行GC含量分析。圖8橫坐標(biāo)是GC含量，縱坐標(biāo)是平均測序深度。以10 kp為窗口無重復(fù)進(jìn)行計(jì)算。根據(jù)這個(gè)圖可以分析這個(gè)物種的GC含量，可以對該樣品是否有外源DNA污染進(jìn)行判斷。
[0048]如圖9所示橫坐標(biāo)為GC含量，縱坐標(biāo)為該GC下的窗口數(shù)目所占比例。將組裝的基因組序列以500 bp為窗口計(jì)算GC含量，窗口之間有250 bp重疊。如圖9示出四種物種(Panda、Dog、Human和Mouse)的基因組GC含量分布分析可以比較組裝基因組與近緣物種的GC含量分布存在的差異。一般情況下，近緣物種的GC分布有相似的趨勢，通過這個(gè)分布圖，可以初步判斷兩個(gè)物種的基因組序列差異。
[0049]在一個(gè)實(shí)施例中，如圖10所示，對于測序深度分析，完成基因組組裝后，為了評價(jià)單堿基覆蓋深度并反映堿基準(zhǔn)確性，進(jìn)行堿基深度分析。如圖10橫坐標(biāo)為深度，縱坐標(biāo)為該深度的堿基數(shù)目占所有堿基數(shù)目的比例。采用比對程序SOAPaligner將過濾之后的Reads比對回拼接的基因組序列上，然后根據(jù)比對結(jié)果統(tǒng)計(jì)每個(gè)堿基被覆蓋的次數(shù)，從而可以得到各種測序深度的堿基占全基因組的百分比。
[0050]在一個(gè)實(shí)施例中，對于常染色體區(qū)域覆蓋度評估，把已經(jīng)公布或者獲得的BAC或Fosmid克隆序列作為參考，將組裝完成的基因組序列比對回參考序列上(blast-ele-5)，檢查組裝的序列對已知序列的覆蓋度(對于來自于同一個(gè)體的數(shù)據(jù)，要求覆蓋度達(dá)到 95%水平)。這個(gè)評價(jià)即可以從宏觀的角度判斷組裝基因組的完整度，還可以從微觀的角度去檢查基因組序列的錯(cuò)誤組裝。
[0051]在一個(gè)實(shí)施例中，對于基因區(qū)覆蓋度評估，把已公布或者獲得的EST或轉(zhuǎn)錄組序列作為query序列映射到組裝完成的基因組序列上，檢查組裝序列對已知序列的覆蓋度(來自于同一物種的EST或轉(zhuǎn)錄組數(shù)據(jù)一般要求覆蓋度達(dá)到98%水平)。
[0052]在一個(gè)實(shí)施例中，對于基因聚類分析，首先基因家族是由來至一個(gè)祖先基因的一組基因組成?；蚣易宓蔫b定，是進(jìn)化分析很重要的一個(gè)方面。通過同源基因的聚類及基因家族的鑒定分析，可以得到單拷貝基因家族和多拷貝基因家族。這些家族在物種之間都是比較保守的，可用于物種間親緣關(guān)系的分析。還可以得到物種特有的基因和家族，它們可能和物種的特異性表型有關(guān)。其次，在進(jìn)行具體操作時(shí)，除當(dāng)前測序物種外，共選取6到9個(gè)物種進(jìn)行后續(xù)分析。首先過濾基因集，即過濾掉無起始、無終止、提前終止、cds序列不是3的倍數(shù)的基因，每一個(gè)基因選取最長的轉(zhuǎn)錄本進(jìn)行分析。對所有的蛋白質(zhì)，用比對程序BlastP做比對，對動(dòng)植物分別用TreeFam軟件和OrthoMCL軟件進(jìn)行基因家族聚類分析。
[0053]注:Unclustered genes:該物種特有的基因的個(gè)數(shù)(包含該物種未聚類的特異基因)；Unique famiIys:該物種特有的基因家族(這些基因家族只包含該物種的基因)。
[0054]如圖11所示，其中人、狗、牛、小鼠和大鼠代表胎盤哺乳動(dòng)物，其中，圖11中的A是大多數(shù)的基因是同源直系基因的示意圖是同源直系基因家族的示意圖。
[0055]在一個(gè)實(shí)施例中，對于系統(tǒng)發(fā)育分析，首先用每單拷貝基因家族構(gòu)建物種發(fā)育樹。不同物種的變異速率和該物種的個(gè)體大小或者世代周期有關(guān)，較大的個(gè)體、較長的世代周期一般有較慢的變異速率。其次，在具體操作時(shí)，系統(tǒng)發(fā)育樹是根據(jù)單拷貝基因家族來構(gòu)建的。根據(jù)基因家族聚類結(jié)果，選取單拷貝基因，利用基于蛋白序列的muscle比對后，轉(zhuǎn)換為核苷酸序列，連接形成超基因，在此基礎(chǔ)上利用PhyML或MrBayes重構(gòu)構(gòu)建系統(tǒng)發(fā)育樹。
[0056]如圖12所示為利用直系同源基因的四重兼并位點(diǎn)構(gòu)建系統(tǒng)發(fā)育樹的示意圖，其中，每個(gè)分支長度代表中性進(jìn)化速率；樹枝上的數(shù)字代表dN/dS。
[0057]在一個(gè)實(shí)施例中，對于物種分歧時(shí)間估算，首先每個(gè)單拷貝基因家族中的四重簡并位點(diǎn)，通常用以估算分子鐘(替換速率)以及物種間的分化時(shí)間，其中中性替換速率一般用每個(gè)位點(diǎn)每年的變異數(shù)來衡量。其次在具體操作時(shí)，需要確定標(biāo)定時(shí)間:根據(jù)系統(tǒng)發(fā)育樹分析涉及的物種，查詢找到任意兩個(gè)物種之間的標(biāo)定時(shí)間。如果能在文獻(xiàn)中找到化石標(biāo)定時(shí)間，貝1J最終物種分歧時(shí)間估算結(jié)果可信；如果只能在timetree.0rg網(wǎng)站找到已有文獻(xiàn)、網(wǎng)站的估算結(jié)果作為標(biāo)定時(shí)間，則最終結(jié)果僅供參考；如果找不到標(biāo)定時(shí)間，則無法進(jìn)行分歧時(shí)間估算分析。將標(biāo)定時(shí)間加到物種系統(tǒng)發(fā)育樹中，據(jù)超基因序列，利用mcmctree或multidivtime軟件，以BRMC方法估算物種分歧時(shí)間。
[0058]如圖13所示2.le_9、l.3e_9、l.8e_9和2.6e_9代表替換速率，單位是每個(gè)位點(diǎn)每年；97、61 (15-176),14 (3_53)表示算出來的分化年代，單位是百萬年。人和狗的分化年代來至 TimeTree database (http://www.timetree.0rg),用來作為校正的時(shí)間。
[0059] 在一個(gè)實(shí)施例中，對于基因組共線性分析，首先共線性片段指同一個(gè)物種內(nèi)部或者兩個(gè)物種之間，由于復(fù)制(基因組復(fù)制、染色體復(fù)制或者大片段復(fù)制)或者物種分化而產(chǎn)生的大片段的同源性現(xiàn)象。在該同源片段內(nèi)部，基因在功能上以及排列順序上都是保守的。共線性片段中的基因在物種進(jìn)化過程中保持了高度的保守性。共線性分析將海量數(shù)字相互信息用圖的形式可視化的表現(xiàn)出來，對兩組測序數(shù)據(jù)通過比對分析，發(fā)現(xiàn)他們之間存在的共線性關(guān)系。在具體操作時(shí)，先確定I到2個(gè)需要進(jìn)行基因組共線性分析的近緣物種。對于動(dòng)物，進(jìn)行兩個(gè)物種之間的Iastz比對。如果選擇了 2個(gè)近緣物種,還可以接著用multiz進(jìn)行多物種的序列比對，如果是3個(gè)物種的基因組比對，還需要做出圖13所示的分化時(shí)間和替換速率的估算。對于植物，首先進(jìn)行兩個(gè)物種之間的BLASTP比對。選取同源成對(相互blastp的最佳比對結(jié)果)，再用MCscan尋找物種間同源blocks,選取比對較好的blocks進(jìn)行作圖分析物種間的共線性如圖14所示。
[0060]如圖15所示，圖中的點(diǎn)和線表示局部的比對結(jié)果。若有連好的染色體，此處可以使用Circos畫圖，如圖15所示。
[0061]在一個(gè)實(shí)施例中，如圖16-20所示，對于全基因組復(fù)制分析，首先大片段復(fù)制是基因組中長度為1-200 kb的低拷貝區(qū)域。通常境況下這些復(fù)制區(qū)域包含內(nèi)含子和外顯子等基因結(jié)構(gòu)，因此，相對于普通的基因組序列，這些復(fù)制區(qū)域并不能靠先驗(yàn)獲取。大部分的已報(bào)道過的大片段復(fù)制區(qū)域都是用實(shí)驗(yàn)分析的方法獲取的。對人類基因組的分析表明，在人的基因組中有大量散布的大片段復(fù)制區(qū)域，以及串聯(lián)復(fù)制區(qū)域。資料顯示大片段復(fù)制現(xiàn)象可能在基因組和基因進(jìn)化方面起著重要的作用。全基因組比對結(jié)果是比較基因組分析中的一個(gè)重要基礎(chǔ)，它一般用于識別基因組中的功能元件。例如，通過基因組的多序列比對結(jié)果得到的多個(gè)遠(yuǎn)緣物種的同源序列一般暗示著這些序列是保守的，具有一定的生物特性。在操作時(shí)，動(dòng)物大片段復(fù)制研究(WGAC)對當(dāng)前測序物種基因組序列進(jìn)行自身Iastz比對分析，在比對之前，先屏蔽重復(fù)序列，再將比對的序列切分為多個(gè)小的文件。最大允許的比對空隙為lOObp。比對之后,過濾出長度大于lkbp, identity大于90%的,且去除了重疊部分的區(qū)域，這些區(qū)域即檢測出的大片段復(fù)制區(qū)域。植物大片段復(fù)制研究(WGAC)用當(dāng)前測序物種的蛋白序列進(jìn)行自身BLASTP比對，選取同源成對，再用MCscan尋找同源性blocks，根據(jù)共線性blocks畫圖如圖16所示。植物全基因組復(fù)制(WGD)首先根據(jù)物種之間的進(jìn)化關(guān)系，選取I到2個(gè)物種。將所有蛋白合在一起，BLASTP比對。物種內(nèi)部，選取同源基因；物種之間，選取直系同源基因。根據(jù)MCscan結(jié)果,對兩兩同源基因/直系同源進(jìn)行muscle比對,轉(zhuǎn)換為核苷酸，提取4D位點(diǎn)，計(jì)算4Dtv值。取每個(gè)共線性片段的4Dtv值，如圖20所示。
[0062]植物基因組自身共線性比對分析結(jié)果如圖18所示。若有連好的染色體，此處可以使用Ciross畫圖，如圖19所示。
[0063]實(shí)施本發(fā)明的實(shí)施例的雜合基因組處理方法，具有以下優(yōu)勢:
測序文庫類型增加，在本發(fā)明中可構(gòu)建170bp、500bp、800bp、2k、5k、10k、20k和40kb共8種WGS測序文庫，承諾加入了華大科技獨(dú)有的40kb大片段文庫用于組裝。
[0064]測序深度，測序深度提高，Illumina Hiseq高通量測序，測序深度達(dá)到200倍以上，雜合基因組解決方案采用全基因組鳥槍法(WGS)進(jìn)行測序。根據(jù)基因組大小和重復(fù)序列的具體特點(diǎn)，制定詳細(xì)的基因組測序策略，整體測序深度在200倍以上，以保證單堿基的準(zhǔn)確性和基因組的完整性。根據(jù)組裝需要，我們要進(jìn)行從170bp，500bp，800bp到2kb，1kb,20kb，40kb梯度插入片斷的文庫的兩端測序。不同插入片段文庫的兩端測序?qū)⒖邕^更多不同大小和相似度的重復(fù)序列。通過對單堿基的高深度重復(fù)測序(大于200倍覆蓋度)，使單喊基的質(zhì)量大大提聞。
[0065]算法增加:我們會(huì)構(gòu)建不同長度的插入片段文庫，并采用雙末端測序以提高組裝的準(zhǔn)確性。Contig構(gòu)建。我們采用SOAP denovo軟件,針對于Illumina Hiseq2000技術(shù)特點(diǎn)的短序列組裝分析軟件(S0AP軟件)，并已成功應(yīng)用于黃瓜、熊貓、白菜、螞蟻、馬鈴薯、谷子、牦牛、牡蠣等上百個(gè)基因組的組裝，積累了豐富的基因組組裝經(jīng)驗(yàn)，能高效準(zhǔn)確地對Hiseq2000測序得到的海量數(shù)據(jù)進(jìn)行拼接組裝。目前Soapdenovo已經(jīng)升級到Soapdenovo2,較之前的版本有多項(xiàng)功能改進(jìn)，組裝效果更好。在Contig構(gòu)建中，Soapdenovo2具有可選的多k-mer的功能,較之前版本的單一 k_mer組裝,效果更好。使用S0APDenovo2進(jìn)行組裝，參數(shù)使用如下:
S0APdenovo-127mer pregraph _s lib_120X.cfg -o tobacco -K 77 -p 32 -R -d I>pregraph.log
S0APdenovo-127mer contig -g tobacco -R > contig.log
S0APdenovo-127mer map _s lib_120X.cfg -g tobacco -f -p 32 _k 45 > map.logS0APdenovo-127mer scaff -g tobacco -N xxx 2> scaff.log對于Scaffold構(gòu)建，Soapdenovo2增加Contig特性識別功能,對雜合、重復(fù)等區(qū)段做特殊處理，連接效率更高。
[0066]以上所述僅為本發(fā)明的優(yōu)選實(shí)施例而已，并不用于限制本發(fā)明，對于本領(lǐng)域的技術(shù)人員來說，本發(fā)明可以有各種更改和變化。凡在本發(fā)明的精神和原則之內(nèi)，所作的任何修改、等同替換、改進(jìn)等，均應(yīng)包含在本發(fā)明的權(quán)利要求范圍之內(nèi)。
【權(quán)利要求】
1.一種雜合基因組處理方法，其特征在于，包括: 采用全基因組鳥槍法對目標(biāo)物種進(jìn)行至少200倍的超高深度測序，以獲得有效的讀長短序列Reads ； 對所述有效的Reads執(zhí)行組裝和構(gòu)建支架序列Scaffold，以獲得帶有冗余序列的基因組圖譜； 對所述帶有冗余序列的基因組圖譜執(zhí)行雜合識別處理，從而去除雜合區(qū)域中冗余的Scaffold并保留雜合區(qū)域信息以獲得全基因組圖譜。
2.根據(jù)權(quán)利要求1所述的雜合基因組處理方法，其特征在于，所述雜合識別處理包括基于k-mer分布圖來識別所述雜合區(qū)域。
3.根據(jù)權(quán)利要求2所述的雜合基因組處理方法，其特征在于，所述基于k-mer分布圖來識別所述雜合區(qū)域包括: 通過k-mer分布圖，識別所述帶有冗余序列的基因組圖譜中的Unique k-mer區(qū)域,其中，所述Unique k_mer表示在基因組中只出現(xiàn)一次的k_mer,k_mer表示定長為指定值k的短序列，k為正整數(shù)； 根據(jù)所述Unique k_mer區(qū)域來確定所述雜合區(qū)域。
4.根據(jù)權(quán)利要求3所述的雜合基因組處理方法，其特征在于，所述根據(jù)所述Uniquek-mer區(qū)域來確定所述雜合區(qū)域包括: 統(tǒng)計(jì)每條Scaffold上的Unique k-mer的數(shù)量N以及那些在已統(tǒng)計(jì)過的一條Scaffold上出現(xiàn)過的Unique k-mer的數(shù)量M，其中N和M均為正整數(shù)； 如果M/N的值達(dá)到預(yù)定值，則將這兩條Scaffold確定為對應(yīng)的雜合區(qū)域。
5.根據(jù)權(quán)利要求1-4任一所述的雜合基因組處理方法，其特征在于，至少基于組裝的長度來去除雜合區(qū)域中冗余的Scaffold。
6.根據(jù)權(quán)利要求1-4任一所述的雜合基因組處理方法，其特征在于，在構(gòu)建支架序列Scaffold之前,還包括構(gòu)建重疊群序列Contig。
7.根據(jù)權(quán)利要求1-4任一所述的雜合基因組處理方法，其特征在于，還包括對所述全基因組圖譜進(jìn)行基因組注釋和/或基因組進(jìn)化分析。
8.根據(jù)權(quán)利要求7所述的雜合基因組處理方法，其特征在于，所述基因組注釋包括R印eat注釋、基因預(yù)測、基因功能注釋和ncRNA注釋；所述基因組進(jìn)化分析包括基因家族鑒定、物種系統(tǒng)發(fā)育樹構(gòu)建、物種分歧時(shí)間估算、基因組共線性分析和全基因組復(fù)制分析。
9.根據(jù)權(quán)利要求1-4任一項(xiàng)所述的雜合基因組處理方法，其特征在于，采用全基因組鳥槍法對目標(biāo)物種進(jìn)行至少200倍的超高深度測序，以獲得有效的讀長短序列Reads包括: 構(gòu)建文庫； 采用全基因組鳥槍法對所構(gòu)建的文庫進(jìn)行至少200倍的超高深度測序，以獲得原始的Reads ； 對所述原始的Reads進(jìn)行過濾，以獲得有效的Reads。
10.根據(jù)權(quán)利要求9所述的雜合基因組處理方法，其特征在于，構(gòu)建文庫包括構(gòu)建小片段文庫和構(gòu)建大片段文庫，其中，小片段文庫包括長度不同的多個(gè)插入片段文庫，大片段文庫包括長度不同的多個(gè)配對Mate Pair文庫。
【文檔編號】C12Q1/68GK104164479SQ201410137420
【公開日】2014年11月26日 申請日期:2014年4月4日 優(yōu)先權(quán)日:2014年4月4日
【發(fā)明者】詹東亮 申請人:深圳華大基因科技服務(wù)有限公司