[发明专利]短序列组装中序列片段的过滤方法及系统

说明书

技术领域

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

背景技术

新测序技术产生的短序列有以下两个特点：第一，序列长度短；第二，数据量大。长序列组装常用的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，对序列的组装和后续分析，尤其是大规模数据的分析，大基因组的组装具有重要的意义。研究有效的序列过滤方法，节约内存，提升计算性能成为一个亟待解决的问题。

发明内容

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

接收测序序列；

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

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

计算所述短串频率阈值；

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

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

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

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

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

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

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

优选地，所述阈值为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开始；