一种基于slaf简化基因组测序数据的二次数据识别方法
【技术领域】
[0001] 本发明属于高通量测序技术领域，具体涉及一种基因组测序数据的重新识别方 法，特别是涉及一种基于全局比对算法的测序数据两端index序列的二次数据识别方法。
【背景技术】
[0002] 凭借可靠的测序数据质量和低廉的成本，高通量测序技术作为一种获取生命核苷 酸编码序列的工具，积累了海量的序列数据，为生物信息学发展提供重要数据支持。采用高 通量测序技术开发分子标记正在成为性价比最高的分子标记开发方式。SLAF-seq是北京百 迈客生物科技有限公司研发的一套简化基因组测序技术，旨在通过限制性内切酶结合电泳 切胶的方法筛选出具有特定长度的基因组片段，并针对筛选出的基因组片段采用高通量测 序方式开发分子标记，其在中国专利CN103088120A中已公开。最新版的SLAF2技术将有效 基因组读长扩展到2 X 100bp，并能够一次开发10万个标签以上，获取全基因组范围内最完 整的变异图像（SNPs、InDels)，以实现重要农艺性状功能基因定位的卓越能力。
[0003] 由于SLAF技术的简化能力，使得单个样品的测序成本低廉，十分适于大样本群体 的测序分析。为了增加同lane可识别的样品数，SLAF采用了双端index的识别策略，即不 同样品测序片段两端加组合特异的一对长Sbp的index序列。下机的不同样品的原始数据 混在一起，因此需要根据测序reads两端的index序列识别和分配数据。一次识别采用的 是illumina公司研发的识别程序，该识别程序严格的过滤条件会导致部分数据损失，由于 测序错误，测序reads两端的index序列可能与其已知序列不完全匹配。一次识别的程序在 处理时，只能容许index上存在1个错配，而错配数多于1的数据则不被识别，这种判断标 准由于过于严格，会造成可观的数据损失，特别是测序数据出现黑tile，大部分测序reads 两端会不可避免引入错配，一次识别将导致大部分可用数据的损失，特别是当测序质量存 在严重问题时。为了提高数据的识别效率，避免损失有价值的数据，需要对未识别的数据进 行二次识别，以提高数据的总识别效率。
【发明内容】
[0004] 本发明提供一种基于SLAF简化基因组测序数据的二次数据识别方法，该技术首 先采用优化index组合的方式，确保上机测序的两两index序列之间有3个碱基以上的差 异，以此确保后期二次识别的准确性。在此基础上利用使用全局比对算法结合项目样品安 排的index先验信息对一次识别后的未识别数据进行重新识别。为了降低数据识别的时间 成本，本发明在具体实施中，采用构建畸变index数据库的方式对数据进行二次识别，畸变 index数据库包含所有以前SLAF测序项目产生的发生测序错误而未识别出index序列和其 对应的原始index编号，由于illumina提供index序列有限，不同时期测序的数据都使用 同样的index，因此根据以往积累的数据构建的数据库可用于识别当次的测序数据，这样就 不用每次都逐条对未识别index进行全局比对，大大节约的时间，同时保证识别准确性。
[0005] 具体地，本发明提供的基于SLAF简化基因组测序数据的二次数据识别方法包括 以下步骤：
[0006] (1)优化上机index组合，确保index组合之间的错配数大于3 ;
[0007] (2)构建畸变index数据库；
[0008] (3)针对双端index的SLAF测序数据，根据测序结果中index与预先加入的index 之间的相似性，对测序数据进行一次识别；
[0009] (4)判断一次识别效率是否低于预期值，当低于预期值时，采用构建的畸变index 数据库对步骤（3)中未识别的数据进行二次识别。
[0010] 其中，步骤⑴所述的优化上机index组合的方法为：统计两两上机index序列之 间的错配数，根据两两index序列间的错配数，进行index聚类，聚到一个类中的两两index 之间的错配数高于3。
[0011] 进一步地，本发明采用全局比对算法统计两两上机index序列之间的错配数，计 算方法如下：
[0015] 全局比对算法构建了一个二维矩阵F，矩阵的横轴和纵轴分别代表一条序列上每 个碱基的位置，Fi j代表序列1上第i个位置和序列2上第j个位置的比对得分；在比对 开始前要对矩阵的第一行&和第一列F Μ进行初始化，通过初始化分配给第一行和第一 列的每个位置一个初始得分，分别记为d*i和d*j，d为一个比对空位的罚分；初始化后 矩阵中每个位置上的得分为对角线上前一格位置的得分F 1 u i加上错配罚分或匹配加分 S (AliB,)，对角线所对应的位置正上方一格的位置Fli, i加空位罚分，对角线所对应的位置左 边一格位置F1 u加空位罚分三者中得分最高的一个，如此依次计算每个位置的得分，得到 两条序列的最优比对模式，根据最优比对模式统计两条序列之间的错配数，即得到两两上 机index序列之间的错配数。
[0016] 上述步骤（1)所述的index聚类，是指：以任一 index序列为一个类的起点，逐条 引入其它index序列，确保新引入的index序列与类中的所有index序列之间有3个以上 错配，该过程一直继续直到类中满足要求的index序列数目不再增加。
[0017] 需要注意的是，这960个index是基于申请人现有的1022个由illumina测序公 司提供的index通过上述方法得到，因此当原始index集合改变时，符合3个以上错配要 求的index数也会发生变化，同样的，如果改变序列之间错配数的筛选条件，比如要求两两 index序列间的错配数在4个以上时，得到的符合筛选条件的index个数也会改变。本发明 实施例中，index聚类经过优化组合后，两两index序列之间有3个以上错配的index -共 960个。两两index序列之间有4个以上错配的index -共476个。
[0018] 步骤（2)所述的畸变index数据库包含所有以前SLAF测序项目产生的发生测序 错误而未识别出的index序列和其对应的原始index编号。所述的"以前SLAF测序项目产 生的发生测序错误而未识别出的index序列"是指以往采用SLAF测序积累的测序错误而未 识别出的index序列。
[0019] 本领域技术人员应当理解，在本发明方法的步骤（4)中，所述的预期值通常指一 次识别率不低于70 %，优选不低于80 %，更优选不低于90 %。
[0020] 本发明的方法中，步骤（4)在进行二次识别时，需要为未识别的index序列分配其 所属的样品编号；然后根据index的分类，从未识别数据中分别将属于每个样品的数据提 取出来。
[0021] 进一步地，为未识别的index序列分配其所属的样品编号的方法为：遍历当次测 序数据产生的未识别index序列，当某条未识别序列在畸变数据库中有匹配且该index对 应的原始index被当次的测序样品使用时，根据与原始index序列的错配数确定来自未识 别index的样品归属，就可以为这条未识别序列分配样品编号。
[0022] 更进一步地，为了保证识别结果的可靠性，每条可以匹配到原始index的未识别 index要满足两个条件：1)未识别序列只与唯一一条原始index序列错配数不大于3,反 之，如果某条未识别数据与两条及以上原始index的错配数都不大于3,则认为仍无法确认 该条index的归属；2)在满足1)的基础上，要求匹配到的原始index亦为当次安排上