一种三代PacBio测序数据的比对方法
【专利摘要】本发明提供一种有效降低重复序列造成的比对错误的三代PacBio测序数据的比对方法。它使用二代的Illumina数据建立k?mer模型，提取unique?kmer，在三代PacBio测序数据的比对中，使用这个unique?kmer来作为比对时使用的种子(seed)，能大大地降低重复序列的影响，提高比对的速度。
【专利说明】
-种H代PacB i O测序数据的比对方法
技术领域
[0001] 本发明设及生物信息技术领域，具体设及DNA序列的比对方法，它使用二代的 Illumina测序数据进行建模提取关键信息，并利用运些关键信息来辅助S代化CBio测序数 据的比对。
【背景技术】
[0002] S代化CBiO的测序数据，单次测序的错误率约为15%，专口支持S代的比对软件 并不多，目前使用最多的软件为W下两款:（l)blasr;(2)dalign。
[0003] 运两款都是非常优秀的S代比对软件，能支持化CBiO的高错误率。由于基因组本 身存在重复序列，它们拥有高度相似的序列。而运些比对软件，会将运些重复序列进行比对 和输出，从而影响后续的生物学分析(比如组装，表达量分析等）。
【发明内容】
[0004] 本发明的目的是解决W上提出的问题，提供一种有效降低重复序列造成的比对错 误的S代化CBio测序数据的比对方法。它使用二代的Illumina数据建立kmer模型，提取 unique-kmer，在S代化cBio测序数据的比对中，使用运个unique-kmer来作为比对时使用 的种子(seed)，能大大地降低重复序列的影响，提高比对的速度。
[0005] 本发明是通过W下技术方案实现的：
[0006] 本发明是一种S代化CBiO测序数据的比对方法，它包括W下步骤：
[0007] (1)使用111皿ina测序数据建立kmer模型，从中提取unique-kmer;
[000引 （2)使用unique-mer作为比对的seed进行候选reads筛选
[0009] (3)对候选reads进行详细比对。
[0010] 作为优化，使用jellyfish软件对二代Illumina测序数据进行k-mer统计，根据k- mer分布图获取二倍主峰W内的k-mer作为unique-kmer,并使用比特文件或GATB开源包，对 所述unique-kmer进行存储。
[0011] 作为优化，对于k《17，使用一个大小为2G的比特文件(*. bit)来存储，而对于4> 17的情况，把unique-kmer存入GATB开源包中的(*. h5)文件中。
[001^ 作为优化，在步骤（2)中，使用步骤（1)的unique-kmer，如果reads之间共有的 unique-kmer计数超过3，就把运些reads筛选出来，作为候选reads。
[0013] 作为优化，所述步骤(3)包括W下步骤：
[0014] a.先对比对上的seed进行聚类，算出最可能的比对范围，方法如下：
[0015] 建立坐标系，横坐标代表readl比对上的位置，纵坐标代表read2上比对上的位置， 每个点代表两条read上共有的seed,将运些seed用斜率为1的直线进行聚类，将聚到最多点 的直线作为比对上的区域；
[0016] b.再将比对范围进行小区域分割，对每一个分割区域，使用LCS算法计算相似度， 再对整体进行打分，方法如下：
[0017]假设将比对范围分为n个区域，相似度大于0.8的区域有b个，运些小区域总体的相 似碱基为C个，则区域相似度为b/n，碱基相似度为c/a，最后只保留运两个值都大于0.7的数 据。
[001引本发明的有益效果如下：
[0019] 1、使用二代Illumina测序数据提取unique-kmer，提高比对的准确率和速度。
[0020]在基因组中，存在许多重复序列，有些短重复序列甚至出现成百上千次，从而会影 响比对的准确度，增加比对的时间。为了提高比对的准确度，降低比对时间，我们提取在 contig中只出现一次的k-mer，作为unique-kmer。因为二代Illumina测序数据的质量非常 高，在测序深度足够随机的情况下（一般情况为~40x)，使用Jel Iyfish软件对二代 Illumina测序数据进行kmer统计，可W得到k-mer的分布图（图1)。将峰值2倍内区域的k- mer作为unique-kmer。对于k< = 17,使用一个大小为2G的比特文件（*.bit文件）来存储，而 对于k> 17的情况，使用GATB (开源框架），把unique-kmer存入文件(*. h5文件）。其中所使用 的二代Illumina测序数据质量较高，Jellyfish软件具有多线程运行，速度快，内存消耗小 的优点，保证了整个方法具有较高的数据处理质量，W及明显的处理速度优势；
[0021 ] 2、使用unique-kmer作为比对的seed进行候选reads筛选，节约比对时间，提高比 对速度。
[0022] 因为unique-kmer在概率和理论上，在单倍体的基因组中，只会出现一次，从而能 避免重复序列造成的影响。另一方面，由于避免了重复序列的影响，找到的候选reads准确 度非常高，节约了很多比对时间，大大提高了比对速度。
[0023] 3、对候选的reads进行详细比对，节约了内存和比对时间，提高比对速度。
[0024] 很多比对软件的比对方法，都使用了最长公共子序列化CS)的算法，直接对整体区 域进行LCS计算，对于大于IOOk的比对区域则非常浪费内存和时间。本方法也是使用运个算 法，但是做了两方面的改进：（1)事先对seed的比对关系进行聚类，算出最优的比对范围； (2)分区域进行比对。从而节约了内存和比对时间，提高比对速度。
【附图说明】
[00巧]图l:kme;r分布图
[0026] 将所有的数据打断成长度为k的片断(称为k-mer)，横坐标为在k-mer的频数，纵坐 标为该频数k-mer的种类，将峰值2倍内区域的k-mer作为unique-kmer。
[0027] 图2:计算出比对范围示意图
[00巧]图上的每个点代表两条read上共有的seed,横坐标代表readl比对上的位置，纵坐 标代表read2比对上的位置，将运些seed用斜率为1的直线进行聚类，选出聚类最多的直线， 将运个区域作为比对上的范围。
[00巧]图3:本发明流程图
【具体实施方式】
[0030] 下面结合附图对本发明的实施例进行进一步详细说明：
[0031] 实施例：
[0032] (1)使用二代111皿ina测序数据建立kmer模型，从中提取unique-kmer
[00削使用jellyfish软件对二代11 Iumina测序数据进行k-mer统计，将所有的数据打断 成长度为k的片断(称为k-mer)，横坐标为在k-mer的频数，纵坐标为该频数k-mer的种类。根 据k-mer分布图获取二倍主峰W内的k-mer作为unique-kmer,对于k《17，使用一个大小为 2G的比特文件（*. bi t)来存储，而对于k> 17的情况，把unique-kmer存入GATB开源包中的 (*.h5)文件中。其中，二代Illumina测序数据是指通过Illumina公司测序仪获得的二代测 序数据。
[0034]根据上述方法，编写如下程序，用来提取unique-kmer，具体操作命令使用说明如 下： 「00北1
[0036]
[0037] 具体案例实施操作如下：
[003引从二代的Illumina测序数据中，筛选大约40X的数据，把它写入一个叫fq. 1st文件 中：
[0039]
[0
[0
[0042] 因为选取k = 17,将结果存入比特文件中:kl7.bit
[0043] (2)使用unique-kmer与S代化cbio测序数据进行比对，筛选候选reads
[0044] 使用运个unique-kmer来作为比对时使用的种子（seed)，如果reads间共有的 unique-kmer超过3时，把它们作为候选reads。其中，S代化cbio测序数据是指通过化cbio 公司测序仪获得的二代测序数据。
[0045] 根据上述方法，编写一个比对程序，来对S代化Cbio测序数据进行比对，具体操作 命令使用说明如下：
[0046]
[0047]
[004引具体案例实施操作如下：
[0049] 使用两个S代化cbio测序的数据文件，分别为readl.fa,read2.fa，另外还有一个 二代111皿ina测序数据提取的unique-kmer文件:kl7.bit,运行W下命令来进行比对：
[(K)加 ]
[0化1] (3)对候选reads进行详细比对。
[0化2] a.先对比对上的seed进行聚类，算出最可能的比对范围，方法如下：
[0化3]建立坐标系，横坐标代表readl比对上的位置，纵坐标代表read2上比对上的位置， 每个点代表两条read上共有的seed,将运些seed用斜率为1的直线进行聚类，将聚到最多点 的直线作为比对上的区域；
[0054] b.再将比对范围进行小区域分割（可W设定分割长度为IOObp)，对每一个分割区 域，使用LCS算法计算相似度，再对整体进行打分，方法如下：
[0055] 假设将比对范围分为n个区域，相似度大于0.8的区域有b个，运些小区域总体的相 似碱基为C个，则区域相似度为b/n，碱基相似度为c/a，最后只保留运两个值都大于0.7的数 据。
[0056] W上所述的仅是本发明的优选实施方式，应当指出，对于本技术领域中的普通技 术人员来说，在不脱离本发明核屯、技术特征的前提下，还可W做出若干改进和润饰，运些改 进和润饰也应视为本发明的保护范围。
【主权项】
1. 一种三代PacBio测序数据的比对方法，其特征在于，它包括以下步骤： (1) 使用二代Illumina测序数据建立kmer模型，并从中提取出unique-kmer; (2) 使用unique-kmer把它作为比对的seed，与三代Pacbio测序数据进行比对，筛选出 候选reads; (3) 对候选reads进行详细比对。2. 根据权利要求书1中所述的三代PacBio测序数据的比对方法，其特征在于，在所述步 骤（1)中，使用je 1 lyf ish软件对二代11 lumina测序数据进行k-mer统计，根据k-mer分布图 获取二倍主峰以内的k-mer作为unique-kmer，并使用比特文件或GATB开源包，对所述 unique-kmer进行存储。3. 根据权利要求书2中所述的三代PacBio测序数据的比对方法，其特征在于，对于 17,使用一个大小为2G的比特文件(*· bit)来存储，而对于k> 17的情况，把unique-kmer存 入GATB开源包中的(*. h5)文件中。4. 根据权利要求书1中所述的三代PacBio测序数据的比对方法，其特征在于，在所述步 骤(2)中，使用步骤（1)的unique-kmer，如果reads之间共有的unique-kmer计数超过3,就把 这些reads筛选出来，作为候选reads。5. 根据权利要求书1中所述的三代PacBio测序数据的比对方法，其特征在于，所述步骤 (3)包括以下步骤： a. 先对比对上的seed进行聚类，算出最可能的比对范围，方法如下： 建立坐标系，横坐标代表read 1比对上的位置，纵坐标代表read2上比对上的位置，每个 点代表两条read上共有的seed，将这些seed用斜率为1的直线进行聚类，将聚到最多点的直 线作为比对上的区域； b. 再将比对范围进行小区域分割，对每一个分割区域，使用LCS算法计算相似度，再对 整体进行打分，方法如下： 假设将比对范围分为η个区域，相似度大于0.8的区域有b个，这些小区域总体的相似碱 基为c个，则区域相似度为b/n，碱基相似度为c/a，最后只保留这两个值都大于0.7的数据。
【文档编号】G06F19/20GK106021997SQ201610329027
【公开日】2016年10月12日
【申请日】2016年5月17日
【发明人】詹东亮, 王军, 王军一, 郝美荣, 何荣军, 俞凯成, 高金龙, 蔡庆乐
【申请人】杭州和壹基因科技有限公司