CN105320848B - 一种基于slaf简化基因组测序数据的二次数据识别方法 - Google Patents

Abstract

本发明提供了一种基于SLAF简化基因组测序数据的二次数据识别方法，优化index组合，使上机测序的两两index序列之间有3个碱基以上的差异；使用全局比对算法结合项目样品的index先验信息对一次识别后的未识别数据进行重新识别，通过构建畸变index数据库对数据进行二次识别，畸变index数据库包含所有以前SLAF测序项目产生的发生测序错误而未识别出index序列和其对应的原始index编号，根据以往积累的数据构建的数据库可用于识别当次的测序数据，不用每次都逐条对未识别index进行全局比对，大大节约时间，在确保识别可靠的基础上，通过增加识别的错配数提高数据的识别效率，减少可用数据损失。

Description

一种基于SLAF简化基因组测序数据的二次数据识别方法

技术领域

本发明属于高通量测序技术领域，具体涉及一种基因组测序数据的重新识别方法，特别是涉及一种基于全局比对算法的测序数据两端index序列的二次数据识别方法。

背景技术

凭借可靠的测序数据质量和低廉的成本，高通量测序技术作为一种获取生命核苷酸编码序列的工具，积累了海量的序列数据，为生物信息学发展提供重要数据支持。采用高通量测序技术开发分子标记正在成为性价比最高的分子标记开发方式。SLAF-seq是北京百迈客生物科技有限公司研发的一套简化基因组测序技术，旨在通过限制性内切酶结合电泳切胶的方法筛选出具有特定长度的基因组片段，并针对筛选出的基因组片段采用高通量测序方式开发分子标记，其在中国专利CN103088120A中已公开。最新版的SLAF2技术将有效基因组读长扩展到2×100bp，并能够一次开发10万个标签以上，获取全基因组范围内最完整的变异图像(SNPs、InDels)，以实现重要农艺性状功能基因定位的卓越能力。

由于SLAF技术的简化能力，使得单个样品的测序成本低廉，十分适于大样本群体的测序分析。为了增加同lane可识别的样品数，SLAF采用了双端index的识别策略，即不同样品测序片段两端加组合特异的一对长8bp的index序列。下机的不同样品的原始数据混在一起，因此需要根据测序reads两端的index序列识别和分配数据。一次识别采用的是illumina公司研发的识别程序，该识别程序严格的过滤条件会导致部分数据损失，由于测序错误，测序reads两端的index序列可能与其已知序列不完全匹配。一次识别的程序在处理时，只能容许index上存在1个错配，而错配数多于1的数据则不被识别，这种判断标准由于过于严格，会造成可观的数据损失，特别是测序数据出现黑tile，大部分测序reads两端会不可避免引入错配，一次识别将导致大部分可用数据的损失，特别是当测序质量存在严重问题时。为了提高数据的识别效率，避免损失有价值的数据，需要对未识别的数据进行二次识别，以提高数据的总识别效率。

发明内容

本发明提供一种基于SLAF简化基因组测序数据的二次数据识别方法，该技术首先采用优化index组合的方式，确保上机测序的两两index序列之间有3个碱基以上的差异，以此确保后期二次识别的准确性。在此基础上利用使用全局比对算法结合项目样品安排的index先验信息对一次识别后的未识别数据进行重新识别。为了降低数据识别的时间成本，本发明在具体实施中，采用构建畸变index数据库的方式对数据进行二次识别，畸变index数据库包含所有以前SLAF测序项目产生的发生测序错误而未识别出index序列和其对应的原始index编号，由于illumina提供index序列有限，不同时期测序的数据都使用同样的index，因此根据以往积累的数据构建的数据库可用于识别当次的测序数据，这样就不用每次都逐条对未识别index进行全局比对，大大节约的时间，同时保证识别准确性。

具体地，本发明提供的基于SLAF简化基因组测序数据的二次数据识别方法包括以下步骤：

(1)优化上机index组合，确保index组合之间的错配数大于3；

(2)构建畸变index数据库；

(3)针对双端index的SLAF测序数据，根据测序结果中index与预先加入的index之间的相似性，对测序数据进行一次识别；

(4)判断一次识别效率是否低于预期值，当低于预期值时，采用构建的畸变index数据库对步骤(3)中未识别的数据进行二次识别。

其中，步骤(1)所述的优化上机index组合的方法为：统计两两上机index序列之间的错配数，根据两两index序列间的错配数，进行index聚类，聚到一个类中的两两index之间的错配数高于3。

进一步地，本发明采用全局比对算法统计两两上机index序列之间的错配数，计算方法如下：

F_0j＝d*j

F_i0＝d*i

F_ij＝max(F_i-1,j-1+S(A_i,B_j),F_i,j-1+d,F_i-1,j+d)

全局比对算法构建了一个二维矩阵F，矩阵的横轴和纵轴分别代表一条序列上每个碱基的位置,Fij代表序列1上第i个位置和序列2上第j个位置的比对得分；在比对开始前要对矩阵的第一行F_0j和第一列F_i0进行初始化，通过初始化分配给第一行和第一列的每个位置一个初始得分,分别记为d*i和d*j，d为一个比对空位的罚分；初始化后矩阵中每个位置上的得分为对角线上前一格位置的得分F_i-1,j-1加上错配罚分或匹配加分S(A_i，B_j）,对角线所对应的位置正上方一格的位置F_i,j-1加空位罚分，对角线所对应的位置左边一格位置F_i-1,j加空位罚分三者中得分最高的一个，如此依次计算每个位置的得分，得到两条序列的最优比对模式，根据最优比对模式统计两条序列之间的错配数，即得到两两上机index序列之间的错配数。

上述步骤(1)所述的index聚类，是指：以任一index序列为一个类的起点，逐条引入其它index序列，确保新引入的index序列与类中的所有index序列之间有3个以上错配，该过程一直继续直到类中满足要求的index序列数目不再增加。

需要注意的是，这960个index是基于申请人现有的1022个由illumina测序公司提供的index通过上述方法得到，因此当原始index集合改变时，符合3个以上错配要求的index数也会发生变化，同样的，如果改变序列之间错配数的筛选条件，比如要求两两index序列间的错配数在4个以上时，得到的符合筛选条件的index个数也会改变。本发明实施例中，index聚类经过优化组合后，两两index序列之间有3个以上错配的index一共960个。两两index序列之间有4个以上错配的index一共476个。

步骤(2)所述的畸变index数据库包含所有以前SLAF测序项目产生的发生测序错误而未识别出的index序列和其对应的原始index编号。所述的“以前SLAF测序项目产生的发生测序错误而未识别出的index序列”是指以往采用SLAF测序积累的测序错误而未识别出的index序列。

本领域技术人员应当理解，在本发明方法的步骤(4)中，所述的预期值通常指一次识别率不低于70％，优选不低于80％，更优选不低于90％。

本发明的方法中，步骤(4)在进行二次识别时，需要为未识别的index序列分配其所属的样品编号；然后根据index的分类，从未识别数据中分别将属于每个样品的数据提取出来。

进一步地，为未识别的index序列分配其所属的样品编号的方法为：遍历当次测序数据产生的未识别index序列，当某条未识别序列在畸变数据库中有匹配且该index对应的原始index被当次的测序样品使用时，根据与原始index序列的错配数确定来自未识别index的样品归属，就可以为这条未识别序列分配样品编号。

更进一步地，为了保证识别结果的可靠性，每条可以匹配到原始index的未识别index要满足两个条件：1)未识别序列只与唯一一条原始index序列错配数不大于3，反之，如果某条未识别数据与两条及以上原始index的错配数都不大于3，则认为仍无法确认该条index的归属；2)在满足1)的基础上，要求匹配到的原始index亦为当次安排上机的index，反之若匹配到的原始index并非本次上机所安排的，则说明匹配结果不可靠，因而弃用。

此处的错配数的计算也遵循全局比对算法，计算公式为：

F_0j＝d*j

F_i0＝d*i

F_ij＝max(F_i-1,j-1+S(A_i，B_j)，F_i,j-1+d，F_i-1,j+d)

在本发明的一个实施例中，上机的index原本是不超过960个已知序列，每对index对应一个样品编号。但是在测序中，每条测序reads对应的index序列上不同位置随机产生测序错误，从而与原序列不再吻合。在进行一次识别时，程序会将未识别测序reads对应的index序列输出，从而建立一个未识别index列表。经过一次识别而未产生出的数据都是与原始序列相差2个及以上碱基，发生错误的碱基可以是16bp长的index任意位置发生，其结果是一批上机数据可以产生6～20M条未识别index序列。而要确定每条未识别index序列的归属，需要将每条未识别序列与每个样品对应的原始index序列做全局比对，然后根据与原始index序列的错配数确定来自未识别index的样品归属。虽然全局比对算法适用于处理这种序列长度一致的两两index比对问题，但是20M条左右的index依次与960个原始index做全局比对，全部做完要花费5天左右的时间，因此十分消耗时间。为了解决这个问题，申请人独创了建立畸变数据库的方式来识别数据。由于优化后的index为固定的960个，所有的未识别index序列都是由这960个index发生各种变化产生的，因此往期测序得到的未识别index序列积累到一定量时可以代表绝大多数的未识别序列类型，那么只需根据以往的未识别index构建一个数据库，记录每条未识别index和其对应的原始index，那么在进行当批的二次识别时只需读取数据库，根据数据库中已经分配好的信息就可以给未识别index分配样品编号，二次识别的时间成本可压缩在1小时内，几乎对后期的数据分析周期不造成影响。

本发明的方法中，上述步骤(4)中的从未识别数据中分别将属于每个样品的数据提取出来是指，经过一次识别，所有未识别测序数据都存放在一个文件中，通过每个index序列对应样品编号，根据每条测序reads第一行注释信息中的index序列将测序reads分配给各个样品。

本发明的方法还包括最后根据一、二次识别出的数据对二次识别效率进行统计。

本发明提供了上述方法在高通量基因组测序中的应用。

本发明的有益效果在于：(1)通过index优化，保证index序列之间的区分度，提高一次和二次识别的可靠性。(2)通过二次识别，可将20％～40％的未识别数据重新识别出来，提高数据利用率。(3)利用以往积累的数据建立一个畸变index数据库，通过数据库中记录的未识别index匹配信息对当次下机的未识别数据直接进行所属样品分配，大大降低时间成本，使二次识别的时间成本控制在1小时之内，保证二次识别流程不会影响后期的分析周期。(4)通过限制错配数，最优匹配原始index个数并结合当次实际安排index的信息保证二次识别结果的可靠性。

附图说明

图1为本发明方法的流程图。

具体实施方式

以下实施例进一步说明本发明的内容，但不应理解为对本发明的限制。在不背离本发明精神和实质的情况下，对本发明方法、步骤或条件所作的修改或替换，均属于本发明的范围。

下列实施例中未注明具体条件的实验方法，通常按照常规条件。实施例中的使用index序列由illumina公司提供的原始index序列集合经步骤1所述的优化方法筛选获得，确保两两序列错配数4个以上，以更严格的筛选条件保证本发明的实施效果。一次识别采用的是illumina公司研发的识别程序。

实施例1第161批SLAF测序数据二次识别统计

包括如下步骤：

1、本批数据上机测序前，根据样品数量，设计了index优化组合，采用全局比对算法统计两两上机index序列之间的错配数，计算方法为：

F_0j＝d*j

F_i0＝d*i

F_ij＝max(F_i-1,j-1+S(A_i,B_j),F_i,j-1+d,F_i-1,j+d)

该组合中的476个index选自于对1022个原始index序列，具体筛选方法为：以原index集的任一index序列为一个类的起点，逐条引入其它index序列，确保新引入的index序列与类中的所有index序列之间有4个以上错配，该过程一直继续直到类中满足要求的index序列数目不再增加，经过优化组合后，确保两两index序列之间有4个以上错配的index一共476个。

表1第161批次测序数据的优化index集合

2、构建畸变index数据库。假设在理想情况下，给每个样品分配的index序列在测序过程中不发生测序错误，则根据测序数据的index为其分配所属样品将不存在数据识别的问题。但是在实际生产中用于标示测序数据归属index序列会出现测序错误，这就导致错误的index序列与对应的原始序列不一致，由这种不一致进而产生给数据分配归属的不确定性问题，从而引出了数据识别的问题。然而当固定了每次用于测序index序列类型后，由这些序列产生的错误序列类型也是固定的。因此通过累积多次测序项目产生的错误index序列，就能构建立一个由错误index序列组成的畸变数据库，包括绝大多数的高频错误index类型，并通过全局比对算法为数据库每条错误index分配最优原始index序列来源。在分配原始序列时要满足的条件包括：未识别序列只与唯一一条原始index序列错配数不大于3，反之，如果某条未识别数据与两条及以上原始index的错配数都不大于3，则认为仍无法确认该条index的归属；其次要求匹配到的原始index亦为当次安排上机的index。该畸变数据库建立后，在进行二次识别时只要将数据库读取一遍就可以直接为index序列存在测序错误的数据分配归属，省去了逐条进行序列比对的时间花费，大大提升识别效率。

3、使用Illumina Hiseq 2500测序的数据下机后使用illumina的识别流程对SLAF测序进行一次识别，识别结束后对一次识别效率进行统计，统计结果如下：

表2一次识别效率统计

数据识别在整个测序数据的质量控制检测中处于非常重要的地位，由附录中的测序数据质量流程图可知，完成测序的数据用于分析前要经过一系列的质量检测过程，这些过程包括数据识别，过滤接头，测序质量统计，样品差异分析，混样评估，对照比对效率评估等多个环节，已从各个方面把控数据质量。其中数据识别位于整个质量控制流程的最前端，这一步的效率将直接影响后面的其他流程，只有经过数据识别，确定样品归属的测序数据才能用于分析，而未能识别的数据只能被丢弃。由于每次的上机安排测序量是一定的，当数据识别效率较低时，将导致可用数据偏少，部分项目可能会因为数据量不够而需要重新加测数据，这将使整个项目的分析周期延后2周至3周，花费大量的时间和人力成本。

该批SLAF数据被安排与重测序项目的数据同批上机，重测序数据采用的index接头为固定的一种序列，这造成了测序起始位置碱基分布不随机，影响了测序质量，导致本批下机的大部分测序reads前段数据在对应位置上未能获得碱基类型，记为N。由于一次识别程序只能容许一个错配，因此当序列中出现一个以上N时，一次识别程序表现乏力，导致了大批数据不能识别，因此造成该批次的测序数据一次识别效率只有70％。

3、由于本实施例设定的一次识别率预期值是86％，当低于预期值时，进行二次识别，第161批SLAF测序数据的一次识别率为70％，远低于预期值，因此进入二次识别流程，该流程再次读入一次识别程序未能识别的index序列，并将这些未识别序列与畸变数据库中记录的错误index序列进行匹配，能够完全匹配上的index序列，其对应的原始index序列就从畸变数据库中获得。一次识别完成后，所有未识别测序数据都存放在一个文件中，每条测序reads第一行注释信息记录了index序列，通过上一步明确的错误index和原始index的对应关系信息，和已知的原始index序列和样品编号的对应关系，程序就可以逐条为测序reads分配所属样品，至此完成整个二次识别过程。经过二次识别后，数据识别效率提升至83.4％，而整个二次识别过程耗时仅36分钟，详见表3：

表3二次识别效率统计

实施例2第214批SLAF测序数据二次识别统计

1、设定优化的index组合方法，该批次使用index方案与实施例1中一致，在此不再罗列。

2、构建构建畸变index数据库参照实施例1的方法。

3、使用Illumina Hiseq 2500测序的数据下机后使用illumina的识别流程对SLAF测序进行一次识别，识别结束后对一次识别效率进行统计，统计结果如下：

表4一次识别效率统计

该批SLAF数据是一个纯SLAF测序项目集合，即同一个批次上机的项目都是SLAF项目，因此不会受到其他类型测序数据影响，测序质量正常。一次识别程序效率85.9％，是一个比较正常的水平。

3、由于本实施例设定的一次识别率预期值是85.9％，当低于预期值时，进行二次识别，第214批SLAF测序数据的一次识别率为85.9％，因此进入二次识别流程，经过二次识别后，识别效率提升至89.1％，而整个二次识别过程耗时仅50分钟，详见表5：

表5二次识别效率统计

以上的实施例仅仅是对本发明的优选实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通工程技术人员对本发明的技术方案做出的各种变型和改进，均应落入本发明的权利要求书确定的保护范围内。

Claims (10)

1.一种基于SLAF简化基因组测序数据的二次数据识别方法，包括以下步骤：

(1)优化上机index组合，确保index组合之间的错配数大于3；

(2)构建畸变index数据库；

(3)针对双端index的SLAF测序数据，根据测序结果中index与预先加入的index之间的相似性，对测序数据进行一次识别；

(4)判断一次识别效率是否低于预期值，当低于预期值时，采用构建的畸变index数据库对步骤(3)中未识别的数据进行二次识别。

2.根据权利要求1所述的方法，其特征在于，步骤(1)所述的优化上机index组合的方法为：统计两两上机index序列之间的错配数，根据两两index序列间的错配数，进行index聚类，聚到一个类中的两两index之间的错配数高于3。

3.根据权利要求2所述的方法，其特征在于，采用全局比对算法统计两两上机index序列之间的错配数，计算方法如下：

F_0j＝d*j

F_i0＝d*i

F_ij＝max(F_i-1,j-1+S(A_i，B_j)，F_i,j-1+d，F_i-1,j+d)

全局比对算法构建了一个二维矩阵F，矩阵的横轴和纵轴分别代表一条序列上每个碱基的位置,Fij代表序列1上第i个位置和序列2上第j个位置的比对得分；在比对开始前要对矩阵的第一行F_0j和第一列F_i0进行初始化，通过初始化分配给第一行和第一列的每个位置一个初始得分,分别记为d*i和d*j，d为一个比对空位的罚分；初始化后矩阵中每个位置上的得分为对角线上前一格位置的得分F_i-1,j-1加上错配罚分或匹配加分S(A_i,B_j),对角线所对应的位置正上方一格的位置F_i,j-1加空位罚分，对角线所对应的位置左边一格位置F_i-1,j加空位罚分三者中得分最高的一个，如此依次计算每个位置的得分，得到两条序列的最优比对模式，根据最优比对模式统计两条序列之间的错配数，即得到两两上机index序列之间的错配数。

4.根据权利要求2所述的方法，其特征在于，index聚类方法为：以任一index序列为一个类的起点，逐条引入其它index序列，确保新引入的index序列与类中的所有index序列之间有大于3个错配，该过程一直继续直到类中满足要求的index序列数目不再增加。

5.根据权利要求1-4任一所述的方法，其特征在于，步骤(2)所述的畸变index数据库包含所有以前SLAF测序项目产生的发生测序错误而未识别出的index序列和其对应的原始index编号。

6.根据权利要求1-4任一所述的方法，其特征在于，步骤(4)在进行二次识别时，需要为未识别的index序列分配其所属的样品编号；然后根据index的分类，从未识别数据中分别将属于每个样品的数据提取出来。

7.根据权利要求6所述的方法，其特征在于，为未识别的index序列分配其所属的样品编号的方法为：遍历当次测序数据产生的未识别index序列，当某条未识别序列在畸变数据库中有匹配且该index对应的原始index被当次的测序样品使用时，根据与原始index序列的错配数确定来自未识别index的样品归属，就可以为这条未识别序列分配样品编号。

8.根据权利要求7所述的方法，其特征在于，每条可以匹配到原始index的未识别index要满足两个条件：1)未识别序列只与唯一一条原始index序列错配数不大于3，反之，如果某条未识别数据与两条及以上原始index的错配数都不大于3，则认为仍无法确认该条index的归属；2)在满足1)的基础上，要求匹配到的原始index亦为当次安排上机的index，反之若匹配到的原始index并非本次上机所安排的，则说明匹配结果不可靠，因而弃用。

9.根据权利要求6所述的方法，其特征在于，所述从未识别数据中分别将属于每个样品的数据提取出来是指，经过一次识别，所有未识别测序数据都存放在一个文件中，通过每个index序列对应样品编号，根据每条测序reads第一行注释信息中的index序列将测序reads分配给各个样品。

10.权利要求1-9任一所述的方法在高通量基因组测序中的应用。