CN103065067B - 短序列组装中序列片段的过滤方法及系统 - Google Patents

Abstract

本发明公开了一种短序列组装中序列片段的过滤方法，包括以下步骤：接收测序序列；分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串；将得到的所述短串的序列值及所述短串的出现频率存储为一个节点；计算所述短串频率阈值；将频率小于阈值的短串过滤。本发明还提供了短序列组装中序列片段的过滤系统。本发明的有益效果在于，过滤了错误的短串，减小了组装拼接的短串集合，减小了组装拼接程序所需内存，提高了组装拼接程序的性能；在进行短串节点存储的同时对短串出现的频率进行了统计，操作简单；误差小。

Description

短序列组装中序列片段的过滤方法及系统

技术领域

本发明涉及基因工程技术领域，尤其涉及一种短序列组装中序列片段的过滤方法及系统。

背景技术

新测序技术产生的短序列有以下两个特点：第一，序列长度短；第二，数据量大。长序列组装常用的phrap等软件均为基于序列间的交叠（overlap）来进行拼接组装，此方法运用于短序列上会存在运算量太大的问题，没有实际的应用价值。新兴的短序列组装受到内存、时间等的限制，目前只在较小的原核生物基因组中成功应用。新一代测序分析存在以下难点：第一，海量序列片段，基因组源序列的长度从十万碱基（如猪痘病毒、大肠杆菌）到十亿碱基（如黄种人、黄瓜、熊猫基因组）大小不等，而复杂环境（如海水、人体大肠等）宏基因组数据甚至会达到上百亿碱基，而对这些样本进行测序其覆盖度需达到30倍到100倍，这使得产生的基因序列片段剧增，如亚洲黄种人的基因数据可达到1TB；第二，短序列，随着测序技术的发展，测序读长呈不断减小的趋势，较第一代测序仪的测序长度显著下降，例如454测序仪可以测到400bp,Sanger测序法的测序长度可达1000bp到1200bp；第三，测序错误，在测序产生序列片段的过程中可能伴随由于荧光强度识别问题带来测序误差，例如有可能一个碱基T可能被测序仪读出为A。这些错误是难以避免的，而且这个范围通常是0.5%到2%之间。这就意味着一个长度为75bp的源序列如果带有1%的错误率，那么将导致有一半(1-(1-1%)75=52.9%)的测序产生序列片段可能有错误碱基。针对其中第二个问题，高通量的数据本身就可以生成大规模的k-mer节点，这些节点将被构造成图来分析，而由于测序错误的引入，将使得k-mer节点的数目增大5倍，例如人类基因组测序数据将会产生大约15G的k-mer；由测序错误产生的k-mer，如果进入计算机进行直接处理，将会消耗巨大的内存，例如人类基因组测序数据如果不进行序列过滤清洗的话，将会消耗大约2T的内存来存储这些k-mer所构造的图；测序数据中的错误序列还会在构造的图里面形成错误链接，Tip型错误，泡型错误，这些错误和源基因组序列本身的重复序列，基因突变点位等搅合在一起，这将使得后续的基因序列分析无法进行。因此，在短序列组装前进行过滤，去除错误的k-mer，对序列的组装和后续分析，尤其是大规模数据的分析，大基因组的组装具有重要的意义。研究有效的序列过滤方法，节约内存，提升计算性能成为一个亟待解决的问题。

发明内容

本发明旨在解决上述现有技术中存在的问题，提出一种短序列组装中序列片段的过滤方法，包括以下步骤：

接收测序序列；

分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串；

将得到的所述短串的序列值及所述短串的频率存储为一个节点；

计算所述短串频率阈值；

将频率小于阈值的短串过滤。

优选地，所述节点采用hashmap存储，其中，哈希键为所述序列值，值为所述节点。

优选地，所述将得到的所述短串的序列值及所述短串的频率存储为一个节点的步骤具体为：

根据当前节点的短串的序列值在已存储的节点中查询是否已存有当前节点；

如果没有查询到当前节点，则添加所述当前节点；

如果查询到当前节点，则更新所述节点的频率。

优选地，所述节点中存储短串和互补短串中序列值较大者或较小者。

优选地，所述阈值为T=θ×Cov_R，θ为分类模型参数，Cov_R为测序仪器设定的序列克隆倍数实际值。

优选地，所述计算所述短串频率阈值中包括以下步骤：以短串出现的频率为横坐标，以出现所述频率的短串的个数为纵坐标，绘制频率统计图。

优选地，所述Cov_R的值为所述频率统计图上第一个波峰所在位置对应的覆盖度。

优选地，所述Cov_R的计算方法步骤为：

a、对所有的短串按照出现频率的个数排序，并把短串的个数按频率的大小升序存入一个数组a中；

b、删除数组a中前面递减的短串个数；

c、用数组a的前j个数据求和来初始化Sum0；

d、每次从数组a中取出第i个短串个数，加到Sum_x里面，同时Sum_x减去第i-j个频率短串的个数，其中i大于j且i从j开始；

e、如果Sum_x-1<Sum_x，回到步骤c，直到Sum_x-1>Sum_x，进入下一步骤；

f、用j除以Sum_x，即得到Cov_R。

本发明还提供了一种短序列组装中序列片段的过滤系统，包括：

接收单元，用于接收测序序列；

序列切割单元，用于分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串；

存储统计单元，将得到的所述短串的序列值及所述短串的频率存储为一个节点；

统计计算单元，用于计算所述短串频率阈值；

过滤单元，用于将频率小于阈值的短串过滤。

优选地，所述存储统计单元包括：

查询模块，用于根据得到的短串的序列值在已存储的节点中查询是否已存有当前节点；

节点添加模块，用于在所述查询模块没有查询到当前节点时，添加当前节点；

频率更新模块，用于在所述查询模块查询到当前节点时，更新所述当前节点的频率。

本发明的有益效果在于，过滤了错误的短串，减小了组装拼接的短串集合，减小了组装拼接程序所需内存，提高了组装拼接程序的性能；在进行短串节点存储的同时对短串出现的频率进行了统计，操作简单；误差小。

附图说明

图1是本发明提供的序列片段的过滤方法的实现流程图。

图2是本发明提供的序列片段的过滤的系统的结构图。

图3是本发明实施例中大肠杆菌的测序数据的短串频率统计图。

图4是本发明实施例中变异模型模拟测序数据的短串频率统计图。

图5是本发明实施例中454测序仪模型模拟测序数据的短串频率统计图。

具体实施方式

为了使本领域的技术人员更好的理解本申请的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整的描述。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

在本发明的实施例中，通过分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串（k-mer），并将得到的各短串的序列值存储，统计得到的各所述短串出现的频率，绘制所述短串的频率统计图，计算所述短串频率阈值，将频率小于阈值的短串过滤。

图1所示为本发明实施例提供的短序列组装中序列片段过滤方法的实现流程，详述如下：

在步骤S101中，接收测序序列；

在步骤S102中，分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串（k-mer）；

在步骤S103中，将得到的所述短串的序列值及所述短串的频率存储为一个节点；

在步骤S104中，计算所述短串频率阈值；

在步骤S105中，将频率小于阈值的短串过滤。

在本发明的实施例中，测序序列的碱基长度为25-75，切割成固定碱基长度为21-31的短串。然而，切割得到的短串的长度小于测序序列的长度，其长度可以根据测序序列的长度和实际情况设定。每个节点存储相应短串的序列值和频率。这里，可采用longlongint类型文件存储所述节点，其存储格式如下：

［seq：64，frequency：16，...］；

其中，seq存储短串的序列值，所述序列值的计算方法是使用2位存储一个核苷酸序列，如A用00表示，G用01表示，C用10表示，T用11表示，顺序编码下去生成一个占64位的整数值，并且，考虑到对于偶数长度的短串，其互补短串可能为它本身，例如短串GATC的互补短串为GATC本身。为了防止这种混淆，短串的长度均为奇数，另外，由于本发明实施例中数据结构的限制，短串的长度取21-31的奇数；frequency用16位存储所述短串出现的次数，即频率，频率的取值范围为［0，2¹⁶］；其后面的位数还可以用来存储其他值，例如，可以存储删除标记closed，以标识所述短串是否被删除；也可以存储使用标记in_use，以标识所述短串是否被使用过，还可以存储其他标识。

上述步骤S103具体为：

步骤1，根据当前节点的短串的序列值在已存储的节点中查询是否已存有当前节点；

步骤2，如果没有查询到当前节点，则添加所述当前节点；

步骤3，如果查询到当前节点，则更新所述当前节点的频率。

本发明在存储各节点的同时，对短串的频率进行了统计。在本发明的实施例中，使用hashmap存储各节点，哈希键为序列值，值为节点。例如序列为AAAAA的短串（其互补序列为TTTTT），其序列值为1111111111，频率初始值为1，将其序列值1111111111作为键在hashmap中查询是否已经存有当前节点，如果没有查询到当前节点，则添加所述当前节点存储到hashmap中，其值为所述短串的序列值1111111111，频率初始值为1；如果查询到当前节点，则对所述当前节点频率进行更新，增加1。完成后，执行步骤2，查找下一个短串，直至完成全部短串的查找。

为了降低存储节点所需的空间，作为本发明的一个优选实施例，只用一个节点存储互补的两个短串，节点的序列值取互补的两个短串中较大的序列值。如果一个短串的序列值小于其互补短串的序列值，则节点存储所述互补短串的序列值，例如上例中序列AAAAA的序列值存的就是其互补短串TTTTT的值；如果一个短串的序列值大于其互补短串的序列值，则节点存储所述短串的序列值。当然，节点的序列值也可以存储互补的两个短串中较小的序列值。

当然，也可以用其他结构对各节点进行存储，例如可以用树结构进行存储，使用hashmap存储各节点在内存和使用上与用树状结构存储近似，但是用hashmap存储各节点在访问和修改速度上都明显优于树结构。

步骤S104计算所述短串频率阈值，在本实施例中频率阈值的计算方法如下：

所述阈值为T=θ×Cov_R，θ为分类模型参数，Cov_R为测序仪器设定的序列克隆倍数的实际值。分类模型参数的范围一般在0-10%，当分类模型参数偏小时，被过滤的短串（k-mer）较少，可能保留了更多的错误k-mer；当分类模型参数偏大时，被过滤的短串（k-mer）较多，可能会勿将正确的k-mer也过滤掉了，对后续序列拼接组装或基因分析造成影响。因此，分类模型参数根据实际计算的内存条件，后续序列拼接所使用算法特点等因素进行选择。

测序仪器设定的序列克隆倍数是一个理论值，在实际测序过程中可以设定为某一固定值，但是，由于测序仪的误差和测序过程中的操作误差，测序仪器设定的序列克隆倍数的实际值与理论值相差较大，因此，要根据测序结果对其重新进行计算。

在本发明的一个实施例中，以短串出现的频率为横坐标，出现所述频率的短串的个数为纵坐标绘制频率统计图。根据上述的频率统计图，所述Cov_R的值为所述频率统计图上第一个波峰所在位置对应的覆盖度。

例如，选取大肠杆菌的测序数据进行k-mer频率统计，所述频率统计图横坐标为短串出现的频率，纵坐标为出现所述频率的短串的个数，结果如图3所示，第一个波峰所对应的点为（62，12.68），从图3可读出Cov_R值为62。

在本发明的另一个实施例中，所述Cov_R的值可按如下步骤进行计算：

a、对所有的短串按照出现频率的个数排序，并把短串的个数按频率的大小升序存入一个数组a中；

b、删除数组a中前面递减的短串个数；

c、用数组a的前j个数据求和来初始化Sum0；

d、每次从数组a中取出第i个短串个数，加到Sum_x里面，同时Sum_x减去第i-j个频率短串的个数，其中i大于j且i从j开始；

e、如果Sum_x-1<Sum_x，回到步骤c，直到Sum_x-1>Sum_x，进入下一步骤；

f、用j除以Sum_x，即得到Cov_R。

通过设定的分类模型参数和计算出的测序仪器设定的序列克隆倍数实际值，可以得到频率阈值，将频率小于阈值的短串过滤。

本领域的普通技术人员可以理解，实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，所述的程序可以在存储于一计算机可读取存储介质中，所述的存储介质可以为ROM/RAM、磁盘、光盘等，所述程序用来执行以下步骤：

1，接收测序序列；

2，分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串（k-mer）；

3，将得到的所述短串的序列值及所述短串的频率存储为一个节点；

4，计算所述短串频率阈值；

5，将频率小于阈值的短串过滤。

图2所示为本发明实施例提供的短序列组装中序列片段过滤的系统的结构，为了便于说明仅示出了与本发明实施例相关的部分。

所述短序列组装中序列片段过滤的系统可以用于短序列组装或基因分析中，其中：

接收单元201，用于接收测序序列。

序列切割单元202，用于分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串，其实现方式如上所述，在此不再一一赘述。

存储统计单元203，用于将得到的所述短串的序列值及所述短串的频率存储为一个节点，其实现方式如上所述，在此不再一一赘述。

统计计算单元204，用于计算所述短串频率阈值。

过滤单元205，用于将频率小于阈值的短串过滤。

其中，所述存储统计单元203包括：

查询模块2031，用于根据得到的短串的序列值在已存储的节点中查询是否已存有当前节点。

节点添加模块2032，用于在所述查询模块没有查询到当前节点时，添加当前节点，其实现方式如上所述，不再一一赘述。

频率统计模块2033，用于在所述查询模块查询到当前节点时，更新所述节点的频率，所述节点频率增加1。

以下结合具体的测序仪器模拟数据对本发明的过滤系统进行误差分析。

首先利用变异模型生成的模拟测序数据进行验证。

变异模型：假设一个短序列中每个位置测序仪出错的可能性相同。

令RefSeq的长度为N，并且RefSeq中重叠（repeats）所占的比例为β,测序仪器的误差设定为α，k为denovo拼接算法中所设定的k-mer的长度。

于是，理论上可以得到正确k-mer的个数为K_positive，错误k-mer的个数为K_negative，计算公式分别为

K_positive=N(1-β)

K_negative=k×Cov_R×N×α

最终错误k-mer的个数和正确k-mer的个数的比例是：

$P_{\mathrm{error}}=\frac{K_{\mathrm{negative}}}{K_{\mathrm{positive}}}=\frac{k\&times;{\mathrm{Cov}}_R\&times;\mathrm{\&alpha;}}{(1-\mathrm{\&beta;})}\&GreaterEqual;k\&times;{\mathrm{Cov}}_R\&times;\mathrm{\&alpha;}$

在变异模型下，当Cov_R=30，k=21，α=1%时，根据上述公式可以得到Perror=6.3，即约有86%的k-mer短串是错误的，也就是说，内存将少存储86%的k-mer，从而程序的计算量减少86%。在一般情况下k-mer的错误率是大于80%的。

下面进行实验验证，利用采用变异模型的ProcessData程序生成一套Cov_R=30，k=21，α=1%的模拟测序数据，将上述数据用本发明的过滤系统进行处理，得到的频率统计图见图4。在这套模拟测序数据中，大约生成了1亿5千万个不同的k-mer，其中大约有1亿3千万个k-mer是错误的k-mer，取定θ=1%，通过计算得到实际的Cov_R值为30，于是得出频率阈值为3，通过本发明的过滤程序将所有出现次数小于等于3次的k-mer被认定为错误的k-mer，错误k-mer的数量大约为1亿2800万。于是计算出用本发明的过滤系统处理模拟测序数据的结果为有85%的错误k-mer(模型理论值为86%）。使用CART的混淆表（confusiontable）来进行误差分析（见表1）。

表1.变异模型模拟测序数据误差分析表

从表1可以看出使用本发明的过滤系统处理变异模型模拟测序数据时，该模型的系统误差是1.3%，实际结果是正确但预测结果是错误的kmer个数为0，也就是说并没有丢失正确的k-mer，保留了有用信息，因此不会对后续的基因分析产生影响。但是，要达到一定的正确率，通常需要设定一个偏小的θ，然而为了过滤更多的错误k-mer，需要一个偏大的θ，因此，分类模型参数θ的选取非常重要。

然后利用454测序仪模型生成的模拟测序数据进行验证。

利用采用454测序仪模型的MetaSim程序生成一套Cov_R=30，k=21，α=1%的模拟测序数据，将上述数据用本发明的过滤系统进行处理，得到的频率统计图见图5。在这套模拟测序数据中，大约生成了1亿8700万个不同的k-mer，其中大约有1亿6700万个k-mer是错误的k-mer，取定θ=1%，通过计算得到实际的Cov_R值为30，于是得出频率阈值为3，通过本发明的过滤程序将所有出现次数小于等于3次的k-mer被认定为错误的k-mer，错误k-mer的数量大约为1亿6500万。于是计算出用本发明的过滤系统处理模拟测序数据的结果为有88%的错误k-mer（模型理论值为89%）。使用CART的混淆表（confusiontable）来进行误差分析（见表2）。

表2.454测序仪模型模拟测序数据误差分析表

从表2可以看出使用本发明的过滤系统处理454测序仪模型模拟测序数据时，该模型的系统误差是0.8%，实际结果是正确但预测结果是错误的kmer个数不为0，也就是说丢失了正确的k-mer，丢失了有用信息，会对后续的基因分析产生影响。此时，可以考虑将θ值进行重新设定，如令θ=0.9%进行重新过滤。

以上所述的本发明实施方式，并不构成对本发明保护范围的限定。任何在本发明的精神和原则之内所作的修改、等同替换和改进等，均应包含在本发明的权利要求保护范围之内。

Claims (8)

1.一种短序列组装中序列片段的过滤方法，其特征在于，所述方法包括以下步骤：

接收测序序列；

分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串；

将得到的所述短串的序列值及所述短串的频率存储为一个节点；

计算所述短串频率阈值；

将频率小于阈值的短串过滤；

所述阈值为T＝θ×Cov_R，θ为分类模型参数，Cov_R为测序仪器设定的序列克隆倍数实际值。

2.根据权利要求1所述的过滤方法，其特征在于，所述节点采用hashmap存储，其中，哈希键为所述序列值，值为所述节点。

3.根据权利要求1所述的过滤方法，其特征在于，所述将得到的所述短串的序列值及所述短串的频率存储为一个节点的步骤具体为：

根据当前节点的短串的序列值在已存储的节点中查询是否已存有当前节点；

如果没有查询到当前节点，则添加所述当前节点；

如果查询到当前节点，则更新所述当前节点的频率。

4.根据权利要求1所述的过滤方法，其特征在于，所述节点中存储短串和互补短串中序列值较大者或较小者。

5.根据权利要求1所述的过滤方法，其特征在于，所述计算所述短串频率阈值的步骤包括以下步骤：以短串出现的频率为横坐标，以出现所述频率的短串的个数为纵坐标，绘制频率统计图。

6.根据权利要求5所述的过滤方法，其特征在于，所述Cov_R的值为所述频率统计图上第一个波峰所在位置对应的覆盖度。

7.一种短序列组装中序列片段的过滤系统，其特征在于，所述系统包括：

接收单元，用于接收测序序列；

序列切割单元，用于分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串；

存储统计单元，将得到的所述短串的序列值及所述短串的频率存储为一个节点；

统计计算单元，用于计算所述短串频率阈值；

过滤单元，用于将频率小于阈值的短串过滤；

所述阈值为T＝θ×Cov_R，θ为分类模型参数，Cov_R为测序仪器设定的序列克隆倍数实际值。

8.根据权利要求7所述的系统，其特征在于，所述存储统计单元包括：

查询模块，用于根据得到的短串的序列值在已存储的节点中查询是否已存有当前节点；

节点添加模块，用于在所述查询模块没有查询到当前节点时，添加当前节点；

频率更新模块，用于在所述查询模块查询到当前节点时，更新所述当前节点的频率。