一種三代PacBio測序數(shù)據(jù)的比對方法
【專利摘要】本發(fā)明提供一種有效降低重復序列造成的比對錯誤的三代PacBio測序數(shù)據(jù)的比對方法。它使用二代的Illumina數(shù)據(jù)建立k?mer模型，提取unique?kmer，在三代PacBio測序數(shù)據(jù)的比對中，使用這個unique?kmer來作為比對時使用的種子(seed)，能大大地降低重復序列的影響，提高比對的速度。
【專利說明】
-種H代PacB i O測序數(shù)據(jù)的比對方法
技術領域
[0001] 本發(fā)明設及生物信息技術領域，具體設及DNA序列的比對方法，它使用二代的 Illumina測序數(shù)據(jù)進行建模提取關鍵信息，并利用運些關鍵信息來輔助S代化CBio測序數(shù) 據(jù)的比對。
【背景技術】
[0002] S代化CBiO的測序數(shù)據(jù)，單次測序的錯誤率約為15%，?？谥С諷代的比對軟件 并不多，目前使用最多的軟件為W下兩款:（l)blasr;(2)dalign。
[0003] 運兩款都是非常優(yōu)秀的S代比對軟件，能支持化CBiO的高錯誤率。由于基因組本 身存在重復序列，它們擁有高度相似的序列。而運些比對軟件，會將運些重復序列進行比對 和輸出，從而影響后續(xù)的生物學分析(比如組裝，表達量分析等）。

【發(fā)明內容】

[0004] 本發(fā)明的目的是解決W上提出的問題，提供一種有效降低重復序列造成的比對錯 誤的S代化CBio測序數(shù)據(jù)的比對方法。它使用二代的Illumina數(shù)據(jù)建立kmer模型，提取 unique-kmer，在S代化cBio測序數(shù)據(jù)的比對中，使用運個unique-kmer來作為比對時使用 的種子(seed)，能大大地降低重復序列的影響，提高比對的速度。
[0005] 本發(fā)明是通過W下技術方案實現(xiàn)的：
[0006] 本發(fā)明是一種S代化CBiO測序數(shù)據(jù)的比對方法，它包括W下步驟：
[0007] (1)使用111皿ina測序數(shù)據(jù)建立kmer模型，從中提取unique-kmer;
[000引 （2)使用unique-mer作為比對的seed進行候選reads篩選
[0009] (3)對候選reads進行詳細比對。
[0010] 作為優(yōu)化，使用jellyfish軟件對二代Illumina測序數(shù)據(jù)進行k-mer統(tǒng)計，根據(jù)k- mer分布圖獲取二倍主峰W內的k-mer作為unique-kmer,并使用比特文件或GATB開源包，對 所述unique-kmer進行存儲。
[0011] 作為優(yōu)化，對于k《17，使用一個大小為2G的比特文件(*. bit)來存儲，而對于4> 17的情況，把unique-kmer存入GATB開源包中的(*. h5)文件中。
[001^ 作為優(yōu)化，在步驟（2)中，使用步驟（1)的unique-kmer，如果reads之間共有的 unique-kmer計數(shù)超過3，就把運些reads篩選出來，作為候選reads。
[0013] 作為優(yōu)化，所述步驟(3)包括W下步驟：
[0014] a.先對比對上的seed進行聚類，算出最可能的比對范圍，方法如下：
[0015] 建立坐標系，橫坐標代表readl比對上的位置，縱坐標代表read2上比對上的位置， 每個點代表兩條read上共有的seed,將運些seed用斜率為1的直線進行聚類，將聚到最多點 的直線作為比對上的區(qū)域；
[0016] b.再將比對范圍進行小區(qū)域分割，對每一個分割區(qū)域，使用LCS算法計算相似度， 再對整體進行打分，方法如下：
[0017]假設將比對范圍分為n個區(qū)域，相似度大于0.8的區(qū)域有b個，運些小區(qū)域總體的相 似堿基為C個，則區(qū)域相似度為b/n，堿基相似度為c/a，最后只保留運兩個值都大于0.7的數(shù) 據(jù)。
[001引本發(fā)明的有益效果如下：
[0019] 1、使用二代Illumina測序數(shù)據(jù)提取unique-kmer，提高比對的準確率和速度。
[0020]在基因組中，存在許多重復序列，有些短重復序列甚至出現(xiàn)成百上千次，從而會影 響比對的準確度，增加比對的時間。為了提高比對的準確度，降低比對時間，我們提取在 contig中只出現(xiàn)一次的k-mer，作為unique-kmer。因為二代Illumina測序數(shù)據(jù)的質量非常 高，在測序深度足夠隨機的情況下（一般情況為~40x)，使用Jel Iyfish軟件對二代 Illumina測序數(shù)據(jù)進行kmer統(tǒng)計，可W得到k-mer的分布圖（圖1)。將峰值2倍內區(qū)域的k- mer作為unique-kmer。對于k< = 17,使用一個大小為2G的比特文件（*.bit文件）來存儲，而 對于k> 17的情況，使用GATB (開源框架），把unique-kmer存入文件(*. h5文件）。其中所使用 的二代Illumina測序數(shù)據(jù)質量較高，Jellyfish軟件具有多線程運行，速度快，內存消耗小 的優(yōu)點，保證了整個方法具有較高的數(shù)據(jù)處理質量，W及明顯的處理速度優(yōu)勢；
[0021 ] 2、使用unique-kmer作為比對的seed進行候選reads篩選，節(jié)約比對時間，提高比 對速度。
[0022] 因為unique-kmer在概率和理論上，在單倍體的基因組中，只會出現(xiàn)一次，從而能 避免重復序列造成的影響。另一方面，由于避免了重復序列的影響，找到的候選reads準確 度非常高，節(jié)約了很多比對時間，大大提高了比對速度。
[0023] 3、對候選的reads進行詳細比對，節(jié)約了內存和比對時間，提高比對速度。
[0024] 很多比對軟件的比對方法，都使用了最長公共子序列化CS)的算法，直接對整體區(qū) 域進行LCS計算，對于大于IOOk的比對區(qū)域則非常浪費內存和時間。本方法也是使用運個算 法，但是做了兩方面的改進：（1)事先對seed的比對關系進行聚類，算出最優(yōu)的比對范圍； (2)分區(qū)域進行比對。從而節(jié)約了內存和比對時間，提高比對速度。
【附圖說明】
[00巧]圖l:kme;r分布圖
[0026] 將所有的數(shù)據(jù)打斷成長度為k的片斷(稱為k-mer)，橫坐標為在k-mer的頻數(shù)，縱坐 標為該頻數(shù)k-mer的種類，將峰值2倍內區(qū)域的k-mer作為unique-kmer。
[0027] 圖2:計算出比對范圍示意圖
[00巧]圖上的每個點代表兩條read上共有的seed,橫坐標代表readl比對上的位置，縱坐 標代表read2比對上的位置，將運些seed用斜率為1的直線進行聚類，選出聚類最多的直線， 將運個區(qū)域作為比對上的范圍。
[00巧]圖3:本發(fā)明流程圖
【具體實施方式】
[0030] 下面結合附圖對本發(fā)明的實施例進行進一步詳細說明：
[0031] 實施例：
[0032] (1)使用二代111皿ina測序數(shù)據(jù)建立kmer模型，從中提取unique-kmer
[00削使用jellyfish軟件對二代11 Iumina測序數(shù)據(jù)進行k-mer統(tǒng)計，將所有的數(shù)據(jù)打斷 成長度為k的片斷(稱為k-mer)，橫坐標為在k-mer的頻數(shù)，縱坐標為該頻數(shù)k-mer的種類。根 據(jù)k-mer分布圖獲取二倍主峰W內的k-mer作為unique-kmer,對于k《17，使用一個大小為 2G的比特文件（*. bi t)來存儲，而對于k> 17的情況，把unique-kmer存入GATB開源包中的 (*.h5)文件中。其中，二代Illumina測序數(shù)據(jù)是指通過Illumina公司測序儀獲得的二代測 序數(shù)據(jù)。
[0034]根據(jù)上述方法，編寫如下程序，用來提取unique-kmer，具體操作命令使用說明如 下： 「00北1
[0036]
[0037] 具體案例實施操作如下：
[003引從二代的Illumina測序數(shù)據(jù)中，篩選大約40X的數(shù)據(jù)，把它寫入一個叫fq. 1st文件 中：
[0039]
[0
[0
[0042] 因為選取k = 17,將結果存入比特文件中:kl7.bit
[0043] (2)使用unique-kmer與S代化cbio測序數(shù)據(jù)進行比對，篩選候選reads
[0044] 使用運個unique-kmer來作為比對時使用的種子（seed)，如果reads間共有的 unique-kmer超過3時，把它們作為候選reads。其中，S代化cbio測序數(shù)據(jù)是指通過化cbio 公司測序儀獲得的二代測序數(shù)據(jù)。
[0045] 根據(jù)上述方法，編寫一個比對程序，來對S代化Cbio測序數(shù)據(jù)進行比對，具體操作 命令使用說明如下：
[0046]
[0047]
[004引具體案例實施操作如下：
[0049] 使用兩個S代化cbio測序的數(shù)據(jù)文件，分別為readl.fa,read2.fa，另外還有一個 二代111皿ina測序數(shù)據(jù)提取的unique-kmer文件:kl7.bit,運行W下命令來進行比對：
[(K)加 ]
[0化1] (3)對候選reads進行詳細比對。
[0化2] a.先對比對上的seed進行聚類，算出最可能的比對范圍，方法如下：
[0化3]建立坐標系，橫坐標代表readl比對上的位置，縱坐標代表read2上比對上的位置， 每個點代表兩條read上共有的seed,將運些seed用斜率為1的直線進行聚類，將聚到最多點 的直線作為比對上的區(qū)域；
[0054] b.再將比對范圍進行小區(qū)域分割（可W設定分割長度為IOObp)，對每一個分割區(qū) 域，使用LCS算法計算相似度，再對整體進行打分，方法如下：
[0055] 假設將比對范圍分為n個區(qū)域，相似度大于0.8的區(qū)域有b個，運些小區(qū)域總體的相 似堿基為C個，則區(qū)域相似度為b/n，堿基相似度為c/a，最后只保留運兩個值都大于0.7的數(shù) 據(jù)。
[0056] W上所述的僅是本發(fā)明的優(yōu)選實施方式，應當指出，對于本技術領域中的普通技 術人員來說，在不脫離本發(fā)明核屯、技術特征的前提下，還可W做出若干改進和潤飾，運些改 進和潤飾也應視為本發(fā)明的保護范圍。
【主權項】
1. 一種三代PacBio測序數(shù)據(jù)的比對方法，其特征在于，它包括以下步驟： (1) 使用二代Illumina測序數(shù)據(jù)建立kmer模型，并從中提取出unique-kmer; (2) 使用unique-kmer把它作為比對的seed，與三代Pacbio測序數(shù)據(jù)進行比對，篩選出 候選reads; (3) 對候選reads進行詳細比對。2. 根據(jù)權利要求書1中所述的三代PacBio測序數(shù)據(jù)的比對方法，其特征在于，在所述步 驟（1)中，使用je 1 lyf ish軟件對二代11 lumina測序數(shù)據(jù)進行k-mer統(tǒng)計，根據(jù)k-mer分布圖 獲取二倍主峰以內的k-mer作為unique-kmer，并使用比特文件或GATB開源包，對所述 unique-kmer進行存儲。3. 根據(jù)權利要求書2中所述的三代PacBio測序數(shù)據(jù)的比對方法，其特征在于，對于 17,使用一個大小為2G的比特文件(*· bit)來存儲，而對于k> 17的情況，把unique-kmer存 入GATB開源包中的(*. h5)文件中。4. 根據(jù)權利要求書1中所述的三代PacBio測序數(shù)據(jù)的比對方法，其特征在于，在所述步 驟(2)中，使用步驟（1)的unique-kmer，如果reads之間共有的unique-kmer計數(shù)超過3,就把 這些reads篩選出來，作為候選reads。5. 根據(jù)權利要求書1中所述的三代PacBio測序數(shù)據(jù)的比對方法，其特征在于，所述步驟 (3)包括以下步驟： a. 先對比對上的seed進行聚類，算出最可能的比對范圍，方法如下： 建立坐標系，橫坐標代表read 1比對上的位置，縱坐標代表read2上比對上的位置，每個 點代表兩條read上共有的seed，將這些seed用斜率為1的直線進行聚類，將聚到最多點的直 線作為比對上的區(qū)域； b. 再將比對范圍進行小區(qū)域分割，對每一個分割區(qū)域，使用LCS算法計算相似度，再對 整體進行打分，方法如下： 假設將比對范圍分為η個區(qū)域，相似度大于0.8的區(qū)域有b個，這些小區(qū)域總體的相似堿 基為c個，則區(qū)域相似度為b/n，堿基相似度為c/a，最后只保留這兩個值都大于0.7的數(shù)據(jù)。
【文檔編號】G06F19/20GK106021997SQ201610329027
【公開日】2016年10月12日
【申請日】2016年5月17日
【發(fā)明人】詹東亮, 王軍, 王軍一, 郝美榮, 何榮軍, 俞凱成, 高金龍, 蔡慶樂
【申請人】杭州和壹基因科技有限公司