专利名称：短序列组装中序列片段的过滤方法及系统的制作方法
技术领域：
本发明涉及基因工程技术领域，尤其涉及一种短序列组装中序列片段的过滤方法及系统。
背景技术：
新测序技术产生的短序列有以下两个特点第一，序列长度短；第二，数据量大。长序列组装常用的Phrap等软件均为基于序列间的交叠(overlap)来进行拼接组装，此方法运用于短序列上会存在运算量太大的问题，没有实际的应用价值。新兴的短序列组装受到内存、时间等的限制，目前只在较小的原核生物基因组中成功应用。新一代测序分析存在以下难点第一，海量序列片段，基因组源序列的长度从十万碱基(如猪痘病毒、大肠杆菌)到十亿碱基(如黄种人、黄瓜、熊猫基因组)大小不等，而复杂环境(如海水、人体大肠等)宏基因组数据甚至会达到上百亿碱基，而对这些样本进行测序其覆盖度需达到30倍到100倍，这使得产生的基因序列片段剧增，如亚洲黄种人的基因数据可达到ITB ;第二，短序列，随着测序技术的发展，测序读长呈不断减小的趋势，较第一代测序仪的测序长度显著下降，例如454测序仪可以测到400bp, Sanger测序法的测序长度可达IOOObp到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，对序列的组装和后续分析，尤其是大规模数据的分析，大基因组的组装具有重要的意义。研究有效的序列过滤方法，节约内存，提升计算性能成为一个亟待解决的问题。
发明内容
本发明旨在解决上述现有技术中存在的问题，提出一种短序列组装中序列片段的过滤方法，包括以下步骤接收测序序列；分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串；将得到的所述短串的序列值及所述短串的频率存储为一个节点；计算所述短串频率阈值；
将频率小于阈值的短串过滤。优选地,所述节点采用hash map存储,其中，哈希键为所述序列值,值为所述节点。优选地，所述将得到的所述短串的序列值及所述短串的频率存储为一个节点的步骤具体为根据当前节点的短串的序列值在已存储的节点中查询是否已存有当前节点；如果没有查询到当前节点，则添加所述当前节点；如果查询到当前节点，则更新所述节点的频率。
优选地,所述节点中存储短串和互补短串中序列值较大者或较小者。优选地，所述阈值为T= 0 XCove, 0为分类模型参数，Covk为测序仪器设定的序列克隆倍数实际值。优选地，所述计算所述短串频率阈值中包括以下步骤以短串出现的频率为横坐标，以出现所述频率的短串的个数为纵坐标，绘制频率统计图。优选地，所述Covk的值为所述频率统计图上第一个波峰所在位置对应的覆盖度。优选地,所述Covk的计算方法步骤为a、对所有的短串按照出现频率的个数排序,并把短串的个数按频率的大小升序存入一个数组a中；b、删除数组a中前面递减的短串个数；C、用数组a的前j个数据求和来初始化SumO ；d、每次从数组a中取出第i个短串个数，加到Sumx里面，同时Sumx减去第i_j个频率短串的个数，其中i大于j且i从j开始；e、如果Sunv^Sumx,回到步骤c,直到SunvPSunv进入下一步骤；f、用 j 除以 Sumx,即得到 Covk。本发明还提供了一种短序列组装中序列片段的过滤系统，包括接收单元，用于接收测序序列；序列切割单元，用于分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串；存储统计单元，将得到的所述短串的序列值及所述短串的频率存储为一个节点；统计计算单元，用于计算所述短串频率阈值；过滤单元，用于将频率小于阈值的短串过滤。 优选地,所述存储统计单元包括查询模块，用于根据得到的短串的序列值在已存储的节点中查询是否已存有当前节点；节点添加模块，用于在所述查询模块没有查询到当前节点时，添加当前节点；频率更新模块，用于在所述查询模块查询到当前节点时，更新所述当前节点的频率。本发明的有益效果在于，过滤了错误的短串，减小了组装拼接的短串集合，减小了组装拼接程序所需内存，提高了组装拼接程序的性能；在进行短串节点存储的同时对短串出现的频率进行了统计，操作简单；误差小。
图1是本发明提供的序列片段的过滤方法的实现流程图。图2是本发明提供的序列片段的过滤的系统的结构图。图3是本发明实施例中大肠杆菌的测序数据的短串频率统计图。图4是本发明实施例中变异模型模拟测序数据的短串频率统计图。图5是本发明实施例中454测序仪模型模拟测序数据的短串频率统计图。
具体实施例方式为了使本领域的技术人员更好的理解本申请的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整的描述。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。在本发明的实施例中，通过分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串(k-mer)，并将得到的各短串的序列值存储，统计得到的各所述短串出现的频率，绘制所述短串的频率统计图，计算所述短串频率阈值，将频率小于阈值的短串过滤。图1所示为本发明实施例提供的短序列组装中序列片段过滤方法的实现流程，详述如下在步骤SlOl中，接收测序序列；在步骤S102中，分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串(k-mer)；在步骤S103中，将得到的所述短串的序列值及所述短串的频率存储为一个节点；在步骤S104中，计算所述短串频率阈值；在步骤S105中，将频率小于阈值的短串过滤。在本发明的实施例中，测序序列的碱基长度为25-75，切割成固定碱基长度为21-31的短串。然而，切割得到的短串的长度小于测序序列的长度，其长度可以根据测序序列的长度和实际情况设定。每个节点存储相应短串的序列值和频率。这里，可采用longIongint类型文件存储所述节点,其存储格式如下[seq :64, frequency :16,...]；其中，seq存储短串的序列值，所述序列值的计算方法是使用2位存储一个核苷酸序列，如A用00表不,G用01表不,C用10表不,T用11表不,顺序编码下去生成一个占64位的整数值，并且，考虑到对于偶数长度的短串，其互补短串可能为它本身，例如短串GATC的互补短串为GATC本身。为了防止这种混淆，短串的长度均为奇数，另外，由于本发明实施例中数据结构的限制，短串的长度取21-31的奇数frequency用16位存储所述短串出现的次数，即频率，频率的取值范围为
;其后面的位数还可以用来存储其他值，例如，可以存储删除标记closed，以标识所述短串是否被删除；也可以存储使用标记in_use，以标识所述短串是否被使用过，还可以存储其他标识。上述步骤S103具体为步骤1，根据当前节点的短串的序列值在已存储的节点中查询是否已存有当前节占.步骤2，如果没有查询到当前节点，则添加所述当前节点；
步骤3，如果查询到当前节点，则更新所述当前节点的频率。本发明在存储各节点的同时，对短串的频率进行了统计。在本发明的实施例中，使用hash map存储各节点，哈希键为序列值，值为节点。例如序列为AAAAA的短串(其互补序列为TTTTT)，其序列值为1111111111，频率初始值为1，将其序列值1111111111作为键在hash map中查询是否已经存有当前节点,如果没有查询到当前节点,则添加所述当前节点存储到hash map中，其值为所述短串的序列值1111111111，频率初始值为I ;如果查询到当前节点，则对所述当前节点频率进行更新，增加I。完成后，执行步骤2，查找下一个短串，直至完成全部短串的查找。为了降低存储节点所需的空间，作为本发明的一个优选实施例，只用一个节点存储互补的两个短串，节点的序列值取互补的两个短串中较大的序列值。如果一个短串的序列值小于其互补短串的序列值，则节点存储所述互补短串的序列值，例如上例中序列AAAAA的序列值存的就是其互补短串TTTTT的值；如果一个短串的序列值大于其互补短串的序列值，则节点存储所述短串的序列值。当然，节点的序列值也可以存储互补的两个短串中较小的序列值。当然，也可以用其他结构对各节点进行存储，例如可以用树结构进行存储，使用hash map存储各节点在内存和使用上与用树状结构存储近似,但是用hash map存储各节点在访问和修改速度上都明显优于树结构。步骤S104计算所述短串频率阈值，在本实施例中频率阈值的计算方法如下所述阈值为T= 0 XCove, 0为分类模型参数，Covk为测序仪器设定的序列克隆倍数的实际值。分类模型参数的范围一般在0-10%，当分类模型参数偏小时，被过滤的短串(k-mer)较少,可能保留了更多的错误k_mer ;当分类模型参数偏大时,被过滤的短串(k-mer)较多,可能会勿将正确的k_mer也过滤掉了，对后续序列拼接组装或基因分析造成影响。因此，分类模型参数根据实际计算的内存条件，后续序列拼接所使用算法特点等因素进行选择。测序仪器设定的序列克隆倍数是一个理论值，在实际测序过程中可以设定为某一固定值，但是，由于测序仪的误差和测序过程中的操作误差，测序仪器设定的序列克隆倍数的实际值与理论值相差较大，因此，要根据测序结果对其重新进行计算。在本发明的一个实施例中，以短串出现的频率为横坐标，出现所述频率的短串的个数为纵坐标绘制频率统计图。根据上述的频率统计图，所述Covk的值为所述频率统计图上第一个波峰所在位置对应的覆盖度。例如，选取大肠杆菌的测序数据进行k-mer频率统计，所述频率统计图横坐标为短串出现的频率，纵坐标为出现所述频率的短串的个数，结果如图3所示，第一个波峰所对应的点为(62，12. 68)，从图3可读出Covk值为62。在本发明的另一个实施例中，所述Covk的值可按如下步骤进行计算a、对所有的短串按照出现频率的个数排序，并把短串的个数按频率的大小升序存入一个数组a中；b、删除数组a中前面递减的短串个数；C、用数组a的前j个数据求和来初始化SumO ；d、每次从数组a中取出第i个短串个数，加到Sumx里面，同时Sumx减去第i_j个频率短串的个数，其中i大于j且i从j开始；e、如果Sunv^Sumx,回到步骤c,直到SunvASumx,进入下一步骤；f、用 j 除以 Sumx,即得到 Covk。通过设定的分类模型参数和计算出的测序仪器设定的序列克隆倍数实际值，可以得到频率阈值，将频率小于阈值的短串过滤。本领域的普通技术人员可以理解，实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，所述的程序可以在存储于一计算机可读取存储介质中，所述的存储介质可以为R0M/RAM、磁盘、光盘等，所述程序用来执行以下步骤I，接收测序序列；2，分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串(k-mer)； 3，将得到的所述短串的序列值及所述短串的频率存储为一个节点；4，计算所述短串频率阈值；5，将频率小于阈值的短串过滤。图2所示为本发明实施例提供的短序列组装中序列片段过滤的系统的结构，为了便于说明仅示出了与本发明实施例相关的部分。所述短序列组装中序列片段过滤的系统可以用于短序列组装或基因分析中，其中接收单元201，用于接收测序序列。序列切割单元202，用于分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串，其实现方式如上所述，在此不再一一赘述。存储统计单元203，用于将得到的所述短串的序列值及所述短串的频率存储为一个节点，其实现方式如上所述，在此不再一一赘述。统计计算单元204，用于计算所述短串频率阈值。过滤单元205，用于将频率小于阈值的短串过滤。其中，所述存储统计单元203包括查询模块2031，用于根据得到的短串的序列值在已存储的节点中查询是否已存有当前节点。节点添加模块2032，用于在所述查询模块没有查询到当前节点时，添加当前节点，其实现方式如上所述，不再一一赘述。频率统计模块2033，用于在所述查询模块查询到当前节点时，更新所述节点的频率，所述节点频率增加I。以下结合具体的测序仪器模拟数据对本发明的过滤系统进行误差分析。首先利用变异模型生成的模拟测序数据进行验证。变异模型假设一个短序列中每个位置测序仪出错的可能性相同。令RefSeq的长度为N，并且RefSeq中重叠(r印eats)所占的比例为￠,测序仪器的误差设定为a，k为de novo拼接算法中所设定的k_mer的长度。于是，理论上可以得到正确k-mer的个数为Kptjsitive，错误k-mer的个数为Knegative,计算公式分别为
权利要求
1.一种短序列组装中序列片段的过滤方法，其特征在于，所述方法包括以下步骤 接收测序序列； 分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串； 将得到的所述短串的序列值及所述短串的频率存储为一个节点； 计算所述短串频率阈值； 将频率小于阈值的短串过滤。
2.根据权利要求1所述的过滤方法，其特征在于，所述节点采用hashmap存储,其中，哈希键为所述序列值，值为所述节点。
3.根据权利要求1所述的过滤方法，其特征在于，所述将得到的所述短串的序列值及所述短串的频率存储为一个节点的步骤具体为 根据当前节点的短串的序列值在已存储的节点中查询是否已存有当前节点； 如果没有查询到当前节点，则添加所述当前节点； 如果查询到当前节点，则更新所述当前节点的频率。
4.根据权利要求1所述的过滤方法，其特征在于，所述节点中存储短串和互补短串中序列值较大者或较小者。
5.根据权利要求1所述的过滤方法，其特征在于，所述阈值为T=0 XCove, 0为分类模型参数，Cove为测序仪器设定的序列克隆倍数实际值。
6.根据权利要求5所述的过滤方法，其特征在于，所述计算所述短串频率阈值的步骤包括以下步骤以短串出现的频率为横坐标，以出现所述频率的短串的个数为纵坐标，绘制频率统计图。
7.根据权利要求6所述的过滤方法，其特征在于，所述Covk的值为所述频率统计图上第一个波峰所在位置对应的覆盖度。
8.根据权利要求5所述的过滤方法，其特征在于，所述Covk的计算方法步骤为 a、对所有的短串按照出现频率的个数排序，并把短串的个数按出现频率的大小升序存入一个数组a中； b、删除数组a中前面递减的短串个数； C、用数组a的前j个数据求和来初始化SumO ； d、每次从数组a中取出第i个短串个数，加到Sumx里面，同时Sumx减去第i_j个频率短串的个数，其中i大于j且i从j开始； e、如果Sunv^Sumx,回到步骤c,直到SunvASunv进入下一步骤； f、用j除以Sumx,即得到Covk。
9.一种短序列组装中序列片段的过滤系统，其特征在于，所述系统包括 接收单元，用于接收测序序列； 序列切割单元，用于分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串； 存储统计单元，将得到的所述短串的序列值及所述短串的频率存储为一个节点； 统计计算单元，用于计算所述短串频率阈值； 过滤单元，用于将频率小于阈值的短串过滤。
10.根据权利要求9所述的系统，其特征在于，所述存储统计单元包括查询模块，用于根据得到的短串的序列值在已存储的节点中查询是否已存有当前节点节点添加模块，用于在所述查询模块没有查询到当前节点时，添加当前节点；频率更新模块，用于在所述查询模块查询到当前节点时，更新所述当前节点的频率。
全文摘要
本发明公开了一种短序列组装中序列片段的过滤方法，包括以下步骤接收测序序列；分别将接收到的测序序列逐个碱基滑动切割得到固定碱基长度的短串；将得到的所述短串的序列值及所述短串的出现频率存储为一个节点；计算所述短串频率阈值；将频率小于阈值的短串过滤。本发明还提供了短序列组装中序列片段的过滤系统。本发明的有益效果在于，过滤了错误的短串，减小了组装拼接的短串集合，减小了组装拼接程序所需内存，提高了组装拼接程序的性能；在进行短串节点存储的同时对短串出现的频率进行了统计，操作简单；误差小。
文档编号G06F19/22GK103065067SQ20121057572
公开日2013年4月24日 申请日期2012年12月26日 优先权日2012年12月26日
发明者孟金涛, 魏彦杰, 曾理, 成杰峰, 冯圣中 申请人:深圳先进技术研究院