一種基于基因組組裝的變異檢測方法和裝置的制造方法
【專利摘要】本發(fā)明公開了一種基于基因組組裝的變異檢測方法和裝置，所述方法包括：獲取來源于梯度測序文庫的測序讀段序列；對測序讀段序列進(jìn)行過濾；將已經(jīng)過濾的讀段序列拼接成平均長度達(dá)到第四預(yù)設(shè)值的長序列；將拼接得到的長序列比對到參考基因組上；和對序列比對結(jié)果進(jìn)行變異檢測，獲取序列變異。本發(fā)明的方法有效解決了變異檢測“暗區(qū)”、長序列插入和復(fù)雜結(jié)構(gòu)性變異的檢測難題。
【專利說明】
一種基于基因組組裝的變異檢測方法和裝置
技術(shù)領(lǐng)域
[0001]本發(fā)明涉及基因組學(xué)及生物信息學(xué)技術(shù)領(lǐng)域，具體涉及基于基因組組裝的變異檢測方法和裝置。
【背景技術(shù)】
[0002]隨著國際人類基因組計(jì)劃的完成以及DNA測序成本的急劇下降，借助DNA測序技術(shù)實(shí)現(xiàn)臨床個(gè)性化醫(yī)療的時(shí)代已經(jīng)臨近。在理想的情況下，期望的DNA測序技術(shù)應(yīng)該是一種能夠?qū)⒓?xì)胞中的DNA—次性連續(xù)、完整地讀取出來的技術(shù)。但事實(shí)上，即便是科技如此發(fā)達(dá)的今天也還是做不到，當(dāng)前最先進(jìn)的高通量測序平臺(tái)，也只能將細(xì)胞中完整的DNA序列預(yù)先切割成一小段一小段，然后再在機(jī)器上進(jìn)行測序，并且當(dāng)前每次讀取的序列長度最大也只能達(dá)到250個(gè)堿基對(英文簡稱:bp)，但目前用得最廣還是10bp和150bp這兩種，而每個(gè)人的基因組總共約有30億個(gè)堿基對。因此，在測序完成之后，實(shí)際上只是得到了一大堆非常零碎的小序列而已。
[0003]在這種情況下，為了對測得的人類基因組進(jìn)行研究，往往需要先借助DNA序列比對軟件(如BWA)，將所有這些小序列片段都比對到人類基因組參考序列上，定出每一個(gè)小序列在基因組上的位置，之后才能進(jìn)行下一步的分析和研究。不過，考慮到人與人之間絕大部分DNA序列是相同的，很多時(shí)候，比如針對不同人群的藥物研發(fā)以及即將到來的個(gè)性化醫(yī)療時(shí)代，都不需要研究所有的30億個(gè)堿基對，而只需要關(guān)注存在差異的部分一一包括序列喊基的插入(insert，INS)、缺失(delete，DEL)、重復(fù)(repeat，REP)和倒位(invers1n，INV)等。這種策略可以很好地抓住關(guān)鍵點(diǎn)，同時(shí)也減少了不必要的勞動(dòng)輸出，在科學(xué)研究中是非常有意義的。因此，如何能夠有效、全面地檢測基因組上的差異，發(fā)現(xiàn)人與人之間的不同點(diǎn)，就成了所有這些研究的一個(gè)重要基礎(chǔ)。
[0004]目前，已有很多研究人員和科研機(jī)構(gòu)提出了許多種不同的方案來解決這樣的問題，并發(fā)布了相關(guān)軟件，如GATK、Platypus、Pindel、Dindel、Breakdancer等基因組序列變異檢測軟件。但是它們都存在一些局限性，就是只能檢測到一些長度很短的序列插入和缺失(一般在50bp以下)以及一些發(fā)生了 100bp以上的序列缺失的變異區(qū)域，而無法很有效地解決50bp~1000bp的序列變異、大長度的序列插入、序列倒位、位置變換(Trans)以及其他更為復(fù)雜的序列結(jié)構(gòu)性變異。之所以存在這種局限，是由于:(1)測序的序列讀長(read)比較短，一個(gè)變異要能夠被檢測出來，最有效的證據(jù)就是測得的序列能夠完整地將一個(gè)變異覆蓋過去，但限于當(dāng)前的測序技術(shù)，測序讀長很難滿足這一需求；(2)人類基因組本身存在著約50%的重復(fù)序列，這些重復(fù)序列會(huì)影響比對軟件將短序列定位到基因組上的準(zhǔn)確性；(3)這些軟件本身的算法存在局限性，它們主要是通過統(tǒng)計(jì)學(xué)的方法來推斷這些序列變異的情況，在推斷過程中要考慮測序文庫的片段長度，如Illumina的測序儀，測序文庫構(gòu)建的序列片段長度一般是300bp~500bp，軟件在發(fā)現(xiàn)序列比對的結(jié)果與測序文庫的序列長度存在差異時(shí)才會(huì)判斷可能存在序列變異，但統(tǒng)計(jì)學(xué)的差異性都是有一定適應(yīng)范圍的，這就導(dǎo)致很多變異會(huì)被漏掉，從而形成變異檢測的“暗區(qū)”，比如長度在50bp-1000bp之間的變異是最容易在這類測序文庫中被漏掉的，這些變異并非不存在了，而是檢測的精度在這一長度范圍上變得很差，因此它也就成了一個(gè)序列變異檢測的“暗區(qū)”。而對于那些大長度的序列插入，由于處在變異中的序列根本就沒能比對上參考基因組，這些軟件就更加無能為力了。
[0005]因此，開發(fā)一種全新且準(zhǔn)確性高，并能全面獲得基因組上所有類型變異的方法和流程成為現(xiàn)今后基因組時(shí)代一個(gè)亟待解決的問題。

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

[0006]本發(fā)明提供一種基于基因組組裝的變異檢測方法和裝置，有效解決變異檢測“暗區(qū)”、長序列插入和復(fù)雜結(jié)構(gòu)性變異的檢測難題。
[0007]根據(jù)本發(fā)明的第一方面，本發(fā)明提供一種基于基因組組裝的變異檢測方法，包括:獲取來源于梯度測序文庫的測序讀段序列，該梯度測序文庫為包括至少3個(gè)梯度的測序文庫的集合，上述測序讀段序列為去除樣本接頭序列并進(jìn)行樣本區(qū)分的讀段序列；對上述測序讀段序列進(jìn)行過濾，除去測序質(zhì)量低于第一預(yù)設(shè)值的堿基個(gè)數(shù)超過整條序列堿基個(gè)數(shù)的第二預(yù)設(shè)值的序列，以及序列中測序結(jié)果不確定的堿基個(gè)數(shù)超過整條序列堿基個(gè)數(shù)的第三預(yù)設(shè)值的序列；按照上述梯度測序文庫從小到大的方式，逐步將各個(gè)文庫中已經(jīng)過濾的讀段序列加進(jìn)拼接過程拼接成平均長度達(dá)到第四預(yù)設(shè)值的長序列；將拼接得到的長序列比對到參考基因組上；對序列比對結(jié)果進(jìn)行變異檢測，獲取序列變異。
[0008]根據(jù)本發(fā)明的第二方面，本發(fā)明提供一種基于基因組組裝的變異檢測裝置，包括:數(shù)據(jù)輸入單元，用于輸入數(shù)據(jù)；數(shù)據(jù)輸出單元，用于輸出數(shù)據(jù)；存儲(chǔ)單元，用于存儲(chǔ)數(shù)據(jù)，其中包括可執(zhí)行的程序；處理器，與數(shù)據(jù)輸入單元、數(shù)據(jù)輸出單元及存儲(chǔ)單元數(shù)據(jù)連接，用于執(zhí)行可執(zhí)行的程序，程序的執(zhí)行包括完成上述基于基因組組裝的變異檢測方法。
[0009]根據(jù)本發(fā)明的方法，通過將短小的讀段序列拼接成長序列完整地將一個(gè)變異覆蓋過去，因此有效解決變異檢測“暗區(qū)”、長序列插入和復(fù)雜結(jié)構(gòu)性變異的檢測難題，高效、快速地實(shí)現(xiàn)序列變異的全面檢測，為解碼個(gè)人基因組和實(shí)現(xiàn)個(gè)性化醫(yī)療打下重要基礎(chǔ)，并填補(bǔ)了現(xiàn)有生物信息學(xué)方法在完整檢測基因組序列變異方面的不足。
【附圖說明】
[0010]本發(fā)明的上述和/或附加的方面和優(yōu)點(diǎn)從結(jié)合下面附圖對實(shí)施方式的描述中將變得明顯和容易理解，其中:
圖1為依據(jù)本發(fā)明的一種實(shí)施方式的基于基因組組裝的變異檢測方法流程圖；
圖2為依據(jù)本發(fā)明的一個(gè)實(shí)施例的梯度測序文庫拼接序列的長度變化圖；
圖3為依據(jù)本發(fā)明的一個(gè)實(shí)施例檢測到的基因組序列變異的數(shù)量分布圖；
圖4為依據(jù)本發(fā)明的一個(gè)實(shí)施例檢測到的基因組序列變異的長度分布圖。
【具體實(shí)施方式】
[0011]依據(jù)本發(fā)明的一種實(shí)施方式，提供一種基于基因組組裝的變異檢測方法，參考圖1，包括如下步驟:
S1:獲取來源于梯度測序文庫的測序讀段序列。
[0012]梯度測序文庫通常是將來自目標(biāo)個(gè)體的基因組樣本經(jīng)過打斷成片段，并根據(jù)所選用的測序方法進(jìn)行相應(yīng)的文庫(library)制備獲得的，可選用的測序方法根據(jù)來自的測序平臺(tái)包括但不限于CG(Complete Genomics)、Illumina/ Solexa、ABI/ SOLiD和Roche 454，依據(jù)所選測序平臺(tái)進(jìn)行單端或雙端測序文庫的制備。本發(fā)明一個(gè)實(shí)施例中，以Illumina/Solexa的新一代測序儀Genome Analyzer作為高通量測序方法實(shí)現(xiàn)梯度測序文庫，該測序技術(shù)是一種基于邊合成邊測序(SBS，Sequencing-By-Synthesis)的新型測序方法，通過利用單分子陣列實(shí)現(xiàn)在小型芯片(Flow Cell)上進(jìn)行橋式PCR反應(yīng)。新的可阻斷技術(shù)可實(shí)現(xiàn)每次只合并一個(gè)堿基，不需要標(biāo)記熒光基團(tuán)，再利用相應(yīng)的激光激發(fā)熒光基團(tuán)捕獲激發(fā)光，從而讀取堿基信息。應(yīng)當(dāng)理解，測序方法并不構(gòu)成對本發(fā)明的限制。
[0013]本發(fā)明中梯度測序文庫是指文庫的平均片段長度呈梯度化分布的多個(gè)文庫的集合，例如至少3個(gè)梯度的測序文庫的集合，比如3至7個(gè)平均片段長度逐漸遞增的文庫，優(yōu)選7個(gè)平均片段長度逐漸遞增的文庫，例如片段平均長度為160-200bp、400-600bp、700-1000bp、2kbp、5kbp、1kbp和20kbp共7個(gè)梯度的測序文庫。一般來講，至少應(yīng)該有3個(gè)梯度的測序文庫，比如160-200bp、400-600bp和2kbp這3個(gè)梯度，最好上述7個(gè)梯度的測序文庫均有。梯度測序文庫的數(shù)量和平均片段長度決定了拼接出的長片段的長度和質(zhì)量情況，如果只有160-200bp、400-600bp和2kbp這3個(gè)梯度的小片段文庫，一般N50值達(dá)到10kbp左右；如果上述7個(gè)梯度的測序文庫均有，N50值可以達(dá)到lOMbp。其中N50定義為:用于衡量序列拼接結(jié)果的完整程度，具體是所有拼接出來的各個(gè)長序列片段，按照其長度從大到小排序，然后從最長的序列開始往下累加，當(dāng)累加到某條序列時(shí)，累加出來的長度就達(dá)到或超過所有這些序列的總長的50%時(shí)，這條臨界序列的長度就稱為N50。上述7個(gè)梯度的測序文庫，可以比較有效地保證最后序列拼接的結(jié)果。需要指出的是，本發(fā)明中N50與平均長度是不同的概念。平均長度定義為拼接出來的序列總長度除以序列的條數(shù)。
[0014]梯度測序文庫構(gòu)建過程中的DNA樣本的提取和純化以及接頭連接等操作，可視具體測序平臺(tái)而定，目前已經(jīng)是成熟的技術(shù)。梯度測序文庫的片段大小控制也有專門的儀器可以實(shí)現(xiàn)，例如Covaris的超聲波打斷儀，通過對具體工作參數(shù)的設(shè)置即可得到特定長度的打斷片段。
[0015]梯度測序文庫的測序深度，一般需要滿足小片段(<2kbp)文庫測得深一些，大片段(彡2kbp)文庫測得淺一些，例如160-200bp的測序文庫的測序深度可以是20X~30X，例如 20 X、22 X、25 X、27 X 或 30 X，優(yōu)選 20 X 或 30 X ;400_600bp 和 700_1000bp 的測序文庫的測序深度在20X以下即可，例如18X、15X、12X、10X或8X，優(yōu)選10X~20X ;大片段測序文庫的測序深度一般控制在10X WTd^^n9X、8X、7X、6X、5X或3X，優(yōu)選5 X~10 X。而且，總測序深度最好不要超過80 X，例如70 X~80 X，不然測序錯(cuò)誤率雖然低，但如果測序深度過高，其影響反而也會(huì)被放大。
[0016]依據(jù)測序平臺(tái)的情況，下機(jī)的測序讀段序列已經(jīng)去除了測序接頭序列，并且已經(jīng)依據(jù)樣本接頭序列進(jìn)行了樣本區(qū)分，需要進(jìn)行去除樣本接頭序列并進(jìn)行樣本區(qū)分，是由于樣本接頭序列是在高通量測序中引入的用來區(qū)別不同樣本的序列，由于第二代高通量測序平臺(tái)的測序通量巨大，一般都是將不同樣本混合在一起上機(jī)測序。測序完成以后，需要依據(jù)樣本接頭序列將不同來源的樣本區(qū)分開來，并需要將樣本接頭序列去除，以利于組裝。來源于同一樣本的讀段序列具有相同的樣本接頭序列，而不同來源的讀段序列的樣本接頭序列不同，這些不同的樣本接頭序列組成樣本接頭序列庫。
[0017]本發(fā)明提供一種去除樣本接頭序列并進(jìn)行樣本區(qū)分的方法，包括:去除樣本接頭序列中測序質(zhì)量低于第五預(yù)設(shè)值的堿基個(gè)數(shù)大于第六預(yù)設(shè)值的讀段序列；順序任意地a)將樣本接頭序列與樣本接頭序列庫中的序列進(jìn)行完全匹配，b)假設(shè)樣本接頭序列降解第七預(yù)設(shè)值個(gè)堿基再與樣本接頭序列庫中序列對應(yīng)部分進(jìn)行完全匹配，c)允許樣本接頭序列有第八預(yù)設(shè)值個(gè)堿基的插入進(jìn)行完全匹配，以及d)允許樣本接頭序列有第九預(yù)設(shè)值個(gè)堿基的缺失進(jìn)行完全匹配；舍棄a) ~d)中均無匹配結(jié)果、a)~d)中其中一個(gè)同時(shí)匹配到兩個(gè)結(jié)果以及僅有且同時(shí)c)和d)匹配出結(jié)果的整條序列；將匹配到樣本接頭序列庫中同一序列的讀段序列作為同一樣本來源的序列而實(shí)現(xiàn)樣本區(qū)分，并去除讀段序列中的樣本接頭序列。
[0018]去除樣本接頭序列中測序質(zhì)量低于第五預(yù)設(shè)值的堿基個(gè)數(shù)大于第六預(yù)設(shè)值的讀段序列這一步驟的實(shí)質(zhì)是以樣本接頭序列的狀況作為測序的質(zhì)控指標(biāo)。其中，此處的測序質(zhì)量是衡量測序平臺(tái)的測序結(jié)果的質(zhì)量狀況的指標(biāo)，依據(jù)具體測序平臺(tái)，如本發(fā)明一個(gè)實(shí)施例中的Illumina/ Solexa測序平臺(tái)，可定義為:，其中，_示測序質(zhì)量，
彥示測序錯(cuò)誤率，由測序平臺(tái)根據(jù)測序過程中多個(gè)不同因素給出，越小測序質(zhì)量就越尚。
[0019]依據(jù)具體測序技術(shù)及測序環(huán)境，對第五預(yù)設(shè)值和第六預(yù)設(shè)值進(jìn)行設(shè)置。本發(fā)明經(jīng)研究確定在下列數(shù)值范圍內(nèi)既能去除不合格的序列，同時(shí)又能保證盡可能多的合格序列得以保留，提高測序讀段序列的利用率。具體地，第五預(yù)設(shè)值可以設(shè)置為4~6，第五預(yù)設(shè)值越大，意味著去除標(biāo)準(zhǔn)越嚴(yán)格，相應(yīng)去除的讀段序列越多，在本發(fā)明一個(gè)具體實(shí)施例中第五預(yù)設(shè)值設(shè)置為5 ;第六預(yù)設(shè)值可以設(shè)置為2~4，第六預(yù)設(shè)值越大，意味著去除標(biāo)準(zhǔn)越寬松，相應(yīng)去除的讀段序列越少，在本發(fā)明一個(gè)具體實(shí)施例中第六預(yù)設(shè)值設(shè)置為3。
[0020]上述完全匹配操作包括4個(gè)步驟，即a)、b)、c)和d)，可以按任意順序進(jìn)行。一般步驟a)最先進(jìn)行，因?yàn)檫@一步驟是最嚴(yán)格的匹配，稱得上是真正意義上的“完全匹配”，而步驟b)、c)和d) —般在步驟a)之后進(jìn)行。本發(fā)明的一個(gè)具體實(shí)施例中，優(yōu)選地按a)、b)、c)和d)的順序進(jìn)行匹配操作。
[0021]之所以進(jìn)行步驟b)的匹配，是因?yàn)榭紤]到一系列實(shí)驗(yàn)過程中，樣本接頭序列可能出現(xiàn)降解情況，在降解不嚴(yán)重的情況下，如果只是按照步驟a)最嚴(yán)格的匹配，可能會(huì)使得一部分合格的讀段序列被排出，降低了讀段序列的利用率。依據(jù)具體測序技術(shù)及測序環(huán)境，對步驟b)中的第七預(yù)設(shè)值進(jìn)行設(shè)置。第七預(yù)設(shè)值越大，意味著匹配標(biāo)準(zhǔn)越寬松。本發(fā)明經(jīng)研究確定第七預(yù)設(shè)值在1~2范圍內(nèi)，即假設(shè)樣本接頭序列降解I或2bp，再與樣本接頭序列庫中序列對應(yīng)部分進(jìn)行完全匹配，既能去除不合格的序列，同時(shí)又能保證盡可能多的合格序列得以保留，提高測序讀段序列的利用率。在本發(fā)明一個(gè)具體實(shí)施例中第七預(yù)設(shè)值設(shè)置為
1
[0022]之所以進(jìn)行步驟c)的匹配，是因?yàn)榭紤]到一系列實(shí)驗(yàn)過程中，樣本接頭序列發(fā)生堿基插入，在堿基插入不嚴(yán)重的情況下，如果只是按照步驟a)最嚴(yán)格的匹配，可能會(huì)使得一部分合格的讀段序列被排出，降低了讀段序列的利用率。依據(jù)具體測序技術(shù)及測序環(huán)境，對步驟c)中的第八預(yù)設(shè)值進(jìn)行設(shè)置。第八預(yù)設(shè)值越大，意味著匹配標(biāo)準(zhǔn)越寬松。本發(fā)明經(jīng)研究確定第八預(yù)設(shè)值在1~2范圍內(nèi)，即允許樣本接頭序列有I或2bp堿基的插入進(jìn)行完全匹配，既能去除不合格的序列，同時(shí)又能保證盡可能多的合格序列得以保留，提高測序讀段序列的利用率。在本發(fā)明一個(gè)具體實(shí)施例中第八預(yù)設(shè)值設(shè)置為1，也就是允許樣本接頭序列僅有I個(gè)堿基的插入，從樣本接頭序列起始端進(jìn)行完全匹配操作，當(dāng)出現(xiàn)某堿基無法匹配時(shí)認(rèn)為該堿基為插入堿基，跳過此堿基后繼續(xù)嚴(yán)格的完全匹配操作。
[0023]之所以進(jìn)行步驟d)的匹配，是因?yàn)榭紤]到一系列實(shí)驗(yàn)過程中，樣本接頭序列發(fā)生堿基缺失，在堿基缺失不嚴(yán)重的情況下，如果只是按照步驟a)最嚴(yán)格的匹配，可能會(huì)使得一部分合格的讀段序列被排出，降低了讀段序列的利用率。依據(jù)具體測序技術(shù)及測序環(huán)境，對步驟d)中的第九預(yù)設(shè)值進(jìn)行設(shè)置。第九預(yù)設(shè)值越大，意味著匹配標(biāo)準(zhǔn)越寬松。本發(fā)明經(jīng)研究確定第九預(yù)設(shè)值在1~2范圍內(nèi)，即允許樣本接頭序列有I或2bp堿基的缺失進(jìn)行完全匹配，既能去除不合格的序列，同時(shí)又能保證盡可能多的合格序列得以保留，提高測序讀段序列的利用率。在本發(fā)明一個(gè)具體實(shí)施例中第九預(yù)設(shè)值設(shè)置為1，也就是允許樣本接頭序列僅有I個(gè)堿基的缺失進(jìn)行完全匹配操作。
[0024]完成上述a)、b)、c)和d)四步匹配操作后，可以優(yōu)選地按照a) >b) >c) >d)的優(yōu)先級(jí)順序確定最終的樣本接頭序列的比對結(jié)果，而對于上述四步匹配操作中均無匹配結(jié)果，其中一個(gè)匹配步驟同時(shí)匹配到兩個(gè)結(jié)果，或僅有步驟c)和d)的匹配結(jié)果，認(rèn)為是無效信息，將相應(yīng)的整條讀段序列去除。將匹配到樣本接頭序列庫中同一序列的讀段序列作為同一樣本來源的序列而實(shí)現(xiàn)樣本區(qū)分，并去除讀段序列中的樣本接頭序列。
[0025]S2:對測序讀段序列進(jìn)行過濾。
[0026]由于下機(jī)的測序讀段序列可能包含一些不合格的序列，這些序列的存在會(huì)影響讀段序列的組裝質(zhì)量，因此需要先進(jìn)行過濾處理。
[0027]過濾處理主要是去除以下兩類不合格的序列:測序質(zhì)量低于第一預(yù)設(shè)值的堿基個(gè)數(shù)超過整條序列堿基個(gè)數(shù)的第二預(yù)設(shè)值的序列，以及序列中測序結(jié)果不確定的堿基個(gè)數(shù)(如測序結(jié)果中的N的個(gè)數(shù))超過整條序列堿基個(gè)數(shù)的第三預(yù)設(shè)值的序列。
[0028]測序質(zhì)量是衡量測序平臺(tái)的測序結(jié)果的質(zhì)量狀況的指標(biāo)，依據(jù)具體測序平臺(tái)，如本發(fā)明一個(gè)實(shí)施例中的Illumina/ Solexa測序平臺(tái)，可定義為:。凡彳，其中，ζ表示測序質(zhì)量，凡_彥示測序錯(cuò)誤率，由測序平臺(tái)根據(jù)測序過程中多個(gè)不同因素給出，PerroM小測序質(zhì)量就越高。
[0029]依據(jù)具體測序技術(shù)及測序環(huán)境，對第一預(yù)設(shè)值、第二預(yù)設(shè)值和第三預(yù)設(shè)值進(jìn)行設(shè)置。本發(fā)明經(jīng)研究確定在下列數(shù)值范圍內(nèi)既能去除不合格的序列，同時(shí)又能保證盡可能多的合格序列得以保留，提高測序讀段序列的利用率。具體地，第一預(yù)設(shè)值可以設(shè)置為4~6，第一預(yù)設(shè)值越大，意味著過濾標(biāo)準(zhǔn)越嚴(yán)格，相應(yīng)去除的讀段序列越多，在本發(fā)明一個(gè)具體實(shí)施例中第一預(yù)設(shè)值設(shè)置為5 ;第二預(yù)設(shè)值可以設(shè)置為40%~60%，第二預(yù)設(shè)值越大，意味著過濾標(biāo)準(zhǔn)越寬松，相應(yīng)去除的讀段序列越少，在本發(fā)明一個(gè)具體實(shí)施例中第二預(yù)設(shè)值設(shè)置為50% ;第三預(yù)設(shè)值可以設(shè)置為5%~15%，第三預(yù)設(shè)值越大，意味著過濾標(biāo)準(zhǔn)越寬松，相應(yīng)去除的讀段序列越少，在本發(fā)明一個(gè)具體實(shí)施例中第三預(yù)設(shè)值設(shè)置為10%。
[0030]S3:將已經(jīng)過濾的讀段序列拼接成平均長度達(dá)到第四預(yù)設(shè)值的長序列。
[0031]本發(fā)明的變異檢測方法，在完成以上質(zhì)控(即步驟S2)之后，并不是與以往的變異檢測軟件一樣，直接將這些測序的讀段序列比對到參考基因組上，而是先組裝成長片段，這是本發(fā)明最關(guān)鍵的構(gòu)思之一。將短的讀段序列拼接成長序列，最主要的目的是為了使盡可能多的變異序列能被完整地包含在拼接出來的長序列中，從而解決短序列在這一方面的不足，如【背景技術(shù)】中介紹的變異檢測“暗區(qū)”、長序列插入和復(fù)雜結(jié)構(gòu)性變異的檢測難題。本發(fā)明一個(gè)優(yōu)選實(shí)施例中具體是按照梯度測序文庫從小到大逐步將各個(gè)文庫加進(jìn)拼接過程而實(shí)現(xiàn)長片段拼接。本發(fā)明一個(gè)實(shí)施例中，采用160-200bp、400-600bp、700-1000bp、2kbp、5kbp、1kbp和20kbp共7個(gè)梯度的測序文庫進(jìn)行拼接，使得拼接后的長序列平均長度在I百萬(Mbp)個(gè)堿基對以上，是當(dāng)前最先進(jìn)的二代測序技術(shù)所測到的序列讀長的4000倍！基本上已可以保證人類基因組上絕大部分的序列變異都能被這樣長的序列所包括。用于讀段序列拼接的技術(shù)目前有多種，本發(fā)明一個(gè)實(shí)施例中采用S0APdenOVO2(網(wǎng)址是 http://soap, genomics, org.cn/soapdenov0.html ；參考 Ruibang Luo 等人,“S0APdenovo2: An Empirically Improved Memory-Efficient Short-Read de NovoAssembler., vGigaScience, I (2012)，18〈do1: 10.1186/2047-217X-l_18>)進(jìn)行讀段序列拼接。
[0032]本發(fā)明中的第四預(yù)設(shè)值與梯度測序文庫的數(shù)量以及大小有關(guān)，梯度測序文庫的數(shù)量越多、大小越大，第四預(yù)設(shè)值能夠達(dá)到的數(shù)值越大，例如本發(fā)明一個(gè)實(shí)施例中7個(gè)梯度測序文庫進(jìn)行拼接的結(jié)果是得到長序列平均長度為I百萬個(gè)堿基對。
[0033]S4:將拼接得到的長序列比對到參考基因組上。
[0034]在本發(fā)明的變異檢測方法中，經(jīng)過序列拼接這一步驟之后，將拼接得到的長序列比對到參考基因組上，參考基因組序列可取于公共數(shù)據(jù)庫(例如美國國家生物技術(shù)信息中心(NCBI，nat1nal center for b1technology informat1n)提供的 HG19)。本發(fā)明中，使用拼接后的長序列還可以解決因參考基因組存在重復(fù)性序列而導(dǎo)致短序列出現(xiàn)比對錯(cuò)誤的問題。由于短序列已變成很長的序列，這使得以前用于進(jìn)行短序列比對的軟件(比如BffA)不再適用于本發(fā)明的情景，而只能使用長序列比對軟件，如LAST等。本發(fā)明一個(gè)實(shí)施例中使用LAST軟件(網(wǎng)址是http://last.cbrc.jp/)將這些拼接后的長序列比對到參考基因組HG19上。
[0035]S5:對序列比對結(jié)果進(jìn)行變異檢測，獲取序列變異。
[0036]在本發(fā)明的變異檢測方法中，確定了長序列比對的結(jié)果之后，采用變異檢測技術(shù)識(shí)別變異區(qū)域。主要檢測原理是:掃描整個(gè)序列比對結(jié)果，獲得每一段比對序列中所有與參考基因組存在差異的序列及其位置，檢測出每一比對結(jié)果中的變異；對分段比對的結(jié)果按照拼接序列的物理位置順序排序，挑選出順序不正常的序列段和位置，得到結(jié)構(gòu)性變異區(qū)域；對每個(gè)檢測出來的變異進(jìn)行局部重比對，重新確定其精確的序列位置。本發(fā)明中的該步驟可以采用目前已有的變異檢測軟件實(shí)現(xiàn)，本發(fā)明一個(gè)實(shí)施例中使用AsmVar (網(wǎng)址是http://www.stb1inf.com/AsmVar/)，檢測出所有可能的序列變異。
[0037]由于短序列已經(jīng)拼接成了長序列，使用變異檢測軟件除了能夠檢測到小序列變異，還特別能夠有效地檢測出50bp~1000bp的序列變異、大長度的序列插入、序列倒位、位置變換以及其他更為復(fù)雜的序列結(jié)構(gòu)性變異情況，并可以利用統(tǒng)計(jì)學(xué)的方法判定出這些變異在所研究的樣本中的基因型。
[0038]需要說明的是，上述步驟S3~S5中使用的軟件S0APdenovo2、LAST和AsmVar僅僅是一個(gè)具體實(shí)施例中使用的軟件，并不構(gòu)成對本發(fā)明的限制，因?yàn)槟壳坝衅渌娲夹g(shù)可以實(shí)現(xiàn)同樣的目的。
[0039]本領(lǐng)域普通技術(shù)人員可以理解，上述實(shí)施方式中各種方法的全部或部分步驟可以通過程序來指令相關(guān)硬件完成，該程序可以存儲(chǔ)于一計(jì)算機(jī)可讀存儲(chǔ)介質(zhì)中，存儲(chǔ)介質(zhì)可以包括:只讀存儲(chǔ)器、隨機(jī)存儲(chǔ)器、磁盤或光盤等。
[0040]依據(jù)本發(fā)明的另一方面還提供一種基于基因組組裝的變異檢測裝置，包括:數(shù)據(jù)輸入單元，用于輸入數(shù)據(jù)；數(shù)據(jù)輸出單元，用于輸出數(shù)據(jù)；存儲(chǔ)單元，用于存儲(chǔ)數(shù)據(jù)，其中包括可執(zhí)行的程序；處理器，與上述數(shù)據(jù)輸入單元、數(shù)據(jù)輸出單元及存儲(chǔ)單元數(shù)據(jù)連接，用于執(zhí)行上述可執(zhí)行的程序，該程序的執(zhí)行包括完成上述實(shí)施方式中各種方法的全部或部分步驟。
[0041]以下結(jié)合具體目標(biāo)個(gè)體對依據(jù)本發(fā)明的具體檢測方法的運(yùn)行結(jié)果進(jìn)行詳細(xì)的描述。下述檢測過程所使用的具體參數(shù)設(shè)置為:
1、實(shí)施例樣本:10個(gè)家系共30個(gè)正常人的血液樣本；
2、建立平均長度為160-200bp、400-600bp、700-1000bp、2kbp、5kbp、10kbp 和 20kbp 共 7個(gè)梯度測序文庫，使用Illumina/ Solexa的新一代測序儀Genome Analyzer進(jìn)行基于邊合成邊測序的新型測序方法的測序，其中各個(gè)文庫的測序深度如下:160-200bp文庫測30X，400-600bp文庫測10X，700-1000bp文庫測10X，2kbp~20kb這四個(gè)文庫每個(gè)測5X，最終每個(gè)樣本平均測70 X。
[0042]3、設(shè)置第一預(yù)設(shè)值為5，第二預(yù)設(shè)值為50%，第三預(yù)設(shè)值為10%，第五預(yù)設(shè)值為5，第六預(yù)設(shè)值為3，第七預(yù)設(shè)值、第八預(yù)設(shè)值和第九預(yù)設(shè)值均為1，進(jìn)行過濾、去除樣本接頭序列和樣本區(qū)分，其中匹配操作按照上述a)、b)、c)和d)的順序進(jìn)行；
4、使用S0APdenovo2序列拼接軟件，將過濾之后的短測序讀段序列拼接為精確的長序列，并將所有長度小于10bp的序列過濾掉；
5、采用長序列比對軟件LAST將進(jìn)行過拼接之后的長序列比對到參考基因組上，精確定位這些序列在基因組上的位置；
6、使用AsmVar軟件對整個(gè)序列比對結(jié)果進(jìn)行掃描，獲得所有可能的變異列表，除此之外，AsmVar還對這些變異的基因型進(jìn)行判定，并利用高斯混合模型計(jì)算每個(gè)變異的質(zhì)量值，該值是變異位點(diǎn)可靠性的重要評斷標(biāo)準(zhǔn)。
[0043]本實(shí)施例的主要結(jié)果如下。需要注意的是，這里為了簡明起見給出的圖示只是實(shí)施例的一部分結(jié)果，并非全部。
[0044]圖2示出了在序列拼接過程中，隨著梯度測序文庫的逐步加入，拼接序列的N50隨之不斷增大。序列拼接過程中，嚴(yán)格按照梯度測序文庫從小到大逐步將各個(gè)文庫加進(jìn)拼接的過程中。由圖2所示的結(jié)果可以看到，N50有不斷累積增長的趨勢和過程，直到最后一個(gè)文庫20kbp加進(jìn)去之后，拼接出來的序列N50達(dá)到了 1Mbp左右。圖2中也可以看出，如果不是梯度文庫，比如只有小片段文庫(<2kbp)，那么N50最大也達(dá)不到lOOkbp。另一方面，本發(fā)明中所用的這7個(gè)梯度文庫，其間隔步長是比較合適的，不至于出現(xiàn)前一步與后一步差距過大，導(dǎo)致前后銜接不起來，反而不能起到明顯增長N50的目的。圖2中，附圖符號(hào)“1006-01”等表不30個(gè)正常人的血液樣本中其中一個(gè)的編號(hào)。
[0045]圖3示出的是本實(shí)施例所檢測到的不同序列變異的數(shù)量分布圖。由圖3所示的結(jié)果可以看出，使用本發(fā)明的方法檢測到的序列變異類型更加全面，除了較容易檢測到的序列缺失(即序列刪除)和序列插入之外，還檢測到序列倒位(INV)、位置變換(Trans)和序列片段替換(BSubstitut1n)——包括等長替換和不等長替換。圖3中，附圖符號(hào)“ 1006-01”等表示30個(gè)正常人的血液樣本中其中一個(gè)的編號(hào)。
[0046]圖4示出的是本實(shí)施例檢測到的基因組序列變異的長度分布圖，圖形從中間左右劃分，左邊是序列缺失(即序列刪除)的長度分布，右邊是序列插入的長度分布。需要注意的是，為了作圖的方便，本實(shí)施例將其他類型的變異，包括序列倒位、位置變換以及更為復(fù)雜的序列結(jié)構(gòu)性變異都區(qū)分為序列缺失或者序列插入，具體方法是，如果變異的序列長度小于比對位置斷點(diǎn)之間的長度則區(qū)分為序列缺失，反之則區(qū)分為序列插入。分布圖中的已知轉(zhuǎn)座子數(shù)據(jù)庫、已知短序列插入缺失數(shù)據(jù)庫、單堿基突變數(shù)據(jù)庫、基因組結(jié)構(gòu)性變異數(shù)據(jù)庫和千人數(shù)據(jù)庫表示的柱狀部分是本方法和當(dāng)前短序列結(jié)合統(tǒng)計(jì)推斷的方法均能檢測到的變異集合；而空白柱狀部分(新發(fā)現(xiàn)變異)是本方法檢測到的全新的變異，從中可以非常明顯地看到，本發(fā)明的方法所能檢測到的序列變異比起當(dāng)前短序列結(jié)合統(tǒng)計(jì)推斷的方法，長度范圍更全面，數(shù)量也更多，有效地解決了當(dāng)前序列變異檢測方面的不足。
[0047]以上所述僅為本發(fā)明的較佳實(shí)施例，應(yīng)當(dāng)理解，這些實(shí)施例僅用以解釋本發(fā)明，并不用于限定本發(fā)明。對于本領(lǐng)域的一般技術(shù)人員，依據(jù)本發(fā)明的思想，可以對上述【具體實(shí)施方式】進(jìn)行變化。
【主權(quán)項(xiàng)】
1.一種基于基因組組裝的變異檢測方法，其特征在于，包括: 獲取來源于梯度測序文庫的測序讀段序列，所述梯度測序文庫為包括至少3個(gè)梯度的測序文庫的集合，所述測序讀段序列為去除樣本接頭序列并進(jìn)行樣本區(qū)分的讀段序列； 對所述測序讀段序列進(jìn)行過濾，除去測序質(zhì)量低于第一預(yù)設(shè)值的堿基個(gè)數(shù)超過整條序列堿基個(gè)數(shù)的第二預(yù)設(shè)值的序列，以及序列中測序結(jié)果不確定的堿基個(gè)數(shù)超過整條序列堿基個(gè)數(shù)的第三預(yù)設(shè)值的序列； 按照所述梯度測序文庫從小到大的方式，逐步將各個(gè)文庫中已經(jīng)過濾的讀段序列加進(jìn)拼接過程拼接成平均長度達(dá)到第四預(yù)設(shè)值的長序列； 將拼接得到的長序列比對到參考基因組上； 對序列比對結(jié)果進(jìn)行變異檢測，獲取序列變異。2.根據(jù)權(quán)利要求1所述的方法，其特征在于，所述梯度測序文庫為包括7個(gè)梯度的測序文庫的集合； 優(yōu)選地，所述梯度測序文庫為片段平均長度為160-200bp、400-600bp、700-1000bp、2kbp、5kbp、1kbp和20kbp共7個(gè)梯度的測序文庫的集合。3.根據(jù)權(quán)利要求2所述的方法，其特征在于，160-200bp文庫的測序深度為20X~30X ；400-600bp和700-1OOObp文庫的測序深度分別為20 X以下，優(yōu)選為10X~20X ;2kbp、5kbp、10kbp和20kbp文庫的測序深度分別為1X以下，優(yōu)選為5X-10X ;總測序深度為80X以下，優(yōu)選為70X~80X。4.根據(jù)權(quán)利要求1-3任一項(xiàng)所述的方法，其特征在于，所述測序質(zhì)量通過如下公式計(jì)算，其中，ζ表示測序質(zhì)量，/Wcj彥;示測序錯(cuò)誤率；所述第一預(yù)設(shè)值為4~6，所述第二預(yù)設(shè)值為40%~60%，所述第三預(yù)設(shè)值為5%~15% ;優(yōu)選地，所述第一預(yù)設(shè)值為5，所述第二預(yù)設(shè)值為50% ;所述第三預(yù)設(shè)值為10%。5.根據(jù)權(quán)利要求1-3任一項(xiàng)所述的方法，其特征在于，所述第四預(yù)設(shè)值為50萬以上，優(yōu)選為100萬以上。6.根據(jù)權(quán)利要求1-3任一項(xiàng)所述的方法，其特征在于，所述對序列比對結(jié)果進(jìn)行變異檢測包括: 掃描整個(gè)序列比對結(jié)果，獲得每一段比對序列中所有與參考基因組存在差異的序列及其位置，檢測出每一比對結(jié)果中的變異； 對分段比對的結(jié)果按照拼接序列的物理位置順序排序，挑選出順序不正常的序列段和位置，得到結(jié)構(gòu)性變異區(qū)域； 對每個(gè)檢測出來的變異進(jìn)行局部重比對，重新確定其精確的序列位置。7.根據(jù)權(quán)利要求1-3任一項(xiàng)所述的方法，其特征在于，所述去除樣本接頭序列并進(jìn)行樣本區(qū)分具體包括: 去除樣本接頭序列中測序質(zhì)量低于第五預(yù)設(shè)值的堿基個(gè)數(shù)大于第六預(yù)設(shè)值的讀段序列； 順序任意地a)將樣本接頭序列與樣本接頭序列庫中的序列進(jìn)行完全匹配，b)假設(shè)樣本接頭序列降解第七預(yù)設(shè)值個(gè)堿基再與樣本接頭序列庫中序列對應(yīng)部分進(jìn)行完全匹配，c)允許樣本接頭序列有第八預(yù)設(shè)值個(gè)堿基的插入進(jìn)行完全匹配，以及d)允許樣本接頭序列有第九預(yù)設(shè)值個(gè)堿基的缺失進(jìn)行完全匹配； 舍棄a)~d)中均無匹配結(jié)果、a) ~d)中其中一個(gè)同時(shí)匹配到兩個(gè)結(jié)果以及僅有且同時(shí)c)和d)匹配出結(jié)果的整條序列； 將匹配到樣本接頭序列庫中同一序列的讀段序列作為同一樣本來源的序列而實(shí)現(xiàn)樣本區(qū)分，并去除讀段序列中的樣本接頭序列。8.根據(jù)權(quán)利要求7所述的方法，其特征在于，按順序a)、b)、c)和d)進(jìn)行所述完全匹配。9.根據(jù)權(quán)利要求7或8所述的方法，其特征在于，所述第五預(yù)設(shè)值為4~6，所述第六預(yù)設(shè)值為2~4，所述第七預(yù)設(shè)值、第八預(yù)設(shè)值和第九預(yù)設(shè)值均為1~2 ;優(yōu)選地，所述第五預(yù)設(shè)值為5，所述第六預(yù)設(shè)值為3，所述第七預(yù)設(shè)值、第八預(yù)設(shè)值和第九預(yù)設(shè)值均為I。10.一種基于基因組組裝的變異檢測裝置，其特征在于，包括: 數(shù)據(jù)輸入單元，用于輸入數(shù)據(jù)； 數(shù)據(jù)輸出單元，用于輸出數(shù)據(jù)； 存儲(chǔ)單元，用于存儲(chǔ)數(shù)據(jù)，其中包括可執(zhí)行的程序； 處理器，與所述數(shù)據(jù)輸入單元、數(shù)據(jù)輸出單元及存儲(chǔ)單元數(shù)據(jù)連接，用于執(zhí)行所述可執(zhí)行的程序，所述程序的執(zhí)行包括完成如權(quán)利要求1-9任意一項(xiàng)所述的方法。
【文檔編號(hào)】G06F19/18GK105989246SQ201510043893
【公開日】2016年10月5日
【申請日】2015年1月28日
【發(fā)明人】黃樹嘉, 劉斯洋, 葉偉健, 饒俊華
【申請人】深圳華大基因研究院