多物种未出现k-mer子序列计算和特征分析方法及系统

摘要

本发明提供了一种多物种未出现k‑mer子序列计算和特征分析方法及系统，该方法包括：获取原始的物种全基因组数据，并对物种全基因组数据进行数据预处理，获得全基因组序列预处理数据；基于所述物种全基因组序列预处理数据进行不同k值长度下的未出现k‑mer子序列的计算；对不同的多个物种的全基因组数据，进行多物种之间，相同k值长度下共同未出现k‑mer子序列的计算；基于未出现k‑mer子序列的数据结果，进行长度变量分析、GC含量和AG含量的比例统计分析和差异显著性分析、进行Motif发现。该方法能够有效处理多物种全基因组数据，准确计算出未出现k‑mer子序列，并且做有效分析，且有大大提高计算效率。

说明书

技术领域

本发明涉及一种对多物种全基因组数据进行不同k值长度下的未出现k-mer子序列(即LAUPs)的计算和特征分析的方法。

背景技术

生物序列的k-mer频次统计是生物信息处理中一个非常基础且重要的问题,k-mer频次统计信息可以用来揭示生物序列中各种子序列的分布规律,它是一种衡量序列相似性的重要工具。因而其在物种识别,宏基因组分类,序列拼接,多序列比对及RNA二级结构预测、CpG岛研究等众多的生物学问题上都有着重要的应用。

从国内外研究来看，在k-mer频次计算层面，王树林等研究了k-长DNA子序列在DNA全序列中出现频数的计数问题,设计并实现了k-长DNA子序列内部计数算法和外部计数算法。王磊等针对目前大多数拼接算法对于重复段的处理采用效率较低的反复迭代算法的特点,提出了基于k-mer子串的重复段分析方法,充分考虑了拼接中可能的分割点,设计与分析了识别重复序列并提高序列一致性的高效算法。Carl Kingsford等设计并开发了Jellyfish软件，运用Hash表来存储数据，同时能多线程运行，速度快，内存消耗小，该软件只能运行在64位的Linux系统下。从当前国内外情况来看，k-mer统计的计算方面往往考虑的是出现的排列，以及出现次数。很少有研究考虑到那些未出现的排列。而从在k-mer计算应用层面，RNA二级结构预测是k-mer计算非常热门的方向。例如Tinoco等提出了最小自由能模型,Zuker等针对该模型使用动态规划的方法来寻找最优结构，并且提出了Mfold算法。而这些算法，都需要一个计算的长度参数，这个参数如果过大，算法的复杂度非常高。

k-mer频次统计以及其应用研究是国际研究的热点，并且取得了很多的成果，但是现有技术中仍然存在以下的不足：

(1)大部分研究都从已有的排列入手，没有考虑那些从来未出现的序列。

(2)很少有对那些从来未出现序列的组分特性、排列特性进行的分析和研究。

(3)在k-mer运用在发现新非编码RNA或研究RNA结构预测时，当选取RNA长度过长，算法复杂度非常高。

以下对说明书中出现的英文缩写的含义进行解释如下：

LAUPs：Lineage-associated Underrepresented Permutations，与谱系有关的多物种未出现k-mer排列子序列；

GC：鸟嘌呤(guanine)和胞嘧啶(cytosine)的含量；

AT：腺嘌呤(adenine)和胸腺嘧啶(thymine)的含量；

AG：嘌呤含量，包括腺嘌呤(adenine)和鸟嘌呤(guanine)；

CT：嘧啶含量，包括胞嘧啶(cytosine)和胸腺嘧啶(thymine)。

发明内容

有鉴于现有技术中存在的上述不足，本发明提供了一种针对多物种全基因组数据进行不同k值长度下的未出现k-mer子序列(LAUPs)的计算和特征分析方法，具体而言，本发明提供了如下的技术方案：

一方面，本发明提供了一种多物种未出现k-mer子序列计算和特征分析的方法，所述方法包括：

步骤1、获取原始的物种全基因组数据，提取序列数据，并获得反向互补链数据，所述序列数据与反向互补链数据构成预处理数据；

步骤2、基于所述预处理数据，进行不同k值长度下的未出现k-mer子序列计算；

步骤3、对多个物种的全基因组数据，进行多物种之间，相同k值长度下共同未出现k-mer子序列的计算；

步骤4、基于所述步骤2中未出现k-mer子序列以及步骤3中共同未出现的k-mer子序列的数据结果，进行特性分析。

优选的，所述步骤1中的预处理数据，将序列数据以及反向互补链数据合并成为一个数据文件。

优选的，所述步骤2进一步包括：

步骤201、在不同k值长度下，计算同一个物种已出现k-mer结果；

步骤202、依据k的长度，枚举出所有可能的基因序列排列；

步骤203、将所述步骤202中获得的基因序列排列中的每一个排列，在步骤201中获得的已出现k-mer结果中进行查找，当一排列在已出现k-mer结果中未找到，则将其记为未出现的k-mer排列。

优选的，所述步骤201进一步包括：

步骤2011、在计算物种已出现k-mer结果中，无论序列排列出现多少次，均只记为出现一次；该方法可以通过例如更改计算中的次数计数参数的方式来实现，从而实现批量计算；

步骤2012、对同一个数据进行多个k值的计算，直至完成同一个物种已出现k-mer结果的计算。这样可以有效地减少对于相同的数据的重复输入，从而节省输入输出的时间。

优选的，所述步骤3进一步包括：

在获取了多个步骤2得到的多个物种的未出现k-mer子序列的基础上，计算多个物种的共同未出现k-mer子序列的计算，具体采用如下方式：

其中，CommonLAUPs代表多个物种共同拥有的未出现k-mer子序列，i代表物种的序号，k代表k-mer长度，代表物种i在k值长度下的未出现k-mer排列，N指物种的总数。

优选的，所述步骤4进一步包括：

步骤401、分析首次出现LAUPs的k值大小，得到首次出现的k值阈值范围，并将该阈值范围作为后续计算中的k值阈值；这样可以有效避免对不必要的k值的计算缩小计算中的复杂度；

步骤402、对所述所述步骤2中计算出的未出现k-mer子序列的GC含量和嘌呤含量进行统计；本领域技术人员可以采用常规的统计技术方式，获得GC和嘌呤的统计数据，例如直接进行计数等，只要能获得其统计特性即可；

步骤403、比较未出现k-mer子序列与出现的k-mer子序列的GC含量与AT含量之间、嘌呤含量和嘧啶之间的差异显著性。获得该显著性差异后，即可依据该些数据进行特性分析。

优选的，所述步骤4还包括：

步骤404、对未出现k-mer子序列进行Motif发现分析，获取该些子序列中可能存在的频繁模式。

优选的，所述步骤402中GC含量和嘌呤含量通过以下方式进行统计：

其中，k表示k值的长度，num()函数代表在k长的序列中某种碱基的计数。

另一方面，本发明还提供了一种多物种未出现k-mer子序列计算和特征分析的系统，所述系统包括：

数据预处理模块，用于获取原始的物种全基因组数据，提取序列数据，并获得反向互补链数据，所述序列数据与反向互补链数据构成预处理数据；

同物种计算模块，用于基于所述预处理数据，进行同一物种不同k值长度下的未出现k-mer子序列计算；

多物种计算模块，用于对多个物种的全基因组数据，进行多物种之间，相同k值长度下共同未出现k-mer子序列的计算；

特性分析模块，用于基于所述步骤2中未出现k-mer子序列以及步骤3中共同未出现的k-mer子序列的数据结果，进行特性分析。

优选的，所述特性分析模块进一步包括：

差异显著性计算单元，用于比较未出现k-mer子序列与出现的k-mer子序列的GC含量与AT含量之间、嘌呤含量和嘧啶之间的差异显著性；

频繁模式计算单元，用于通过Motif发现方式，获取未出现k-mer子序列中可能存在的频繁模式。

与现有技术相比，本发明的技术方案将改进后的Jellyfish k-mer计算方法、差异显著性分析方法、Motif发现方法等结合在一起，有效处理多物种全基因组数据，能准确计算出未出现k-mer子序列(LAUPs)，并且做有效分析，而且在计算效率上也有大大的提高。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1为本发明实施例的方法流程图；

图2为本发明实施例的不同k值长度下的未出现k-mer子序列(LAUPs)的计算流程图；

图3为本发明实施例未出现k-mer子序列(LAUPs)的分析方法流程图；

图4为本发明实施例未出现k-mer子序列(LAUPs)差异显著性检验的分析方法流程图。

具体实施方式

下面结合附图对本发明实施例一种应用程序推荐方法及装置进行详细描述。应当明确，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

本领域技术人员应当知晓，下述具体实施例或具体实施方式，是本发明为进一步解释具体的发明内容而列举的一系列优化的设置方式，而该些设置方式之间均是可以相互结合或者相互关联使用的，除非在本发明明确提出了其中某些或某一具体实施例或实施方式无法与其他的实施例或实施方式进行关联设置或共同使用。同时，下述的具体实施例或实施方式仅作为最优化的设置方式，而不作为限定本发明的保护范围的理解。

实施例1：

在一个具体的实施例中，本发明的具体技术方案可通过如下的实施例进行详尽说明。本发明的技术方案包括以下步骤：

一、数据预处理

首先，获取原始的物种全基因组数据，由于全基因组数据通常带有丰富的标记和注释，所以需要从中提取出序列数据。在具体的实施例中，从NCBI数据库中获取物种的.gbff格式的全基因组数据，并且使用Perl语言进行编程获取其中的序列数据，具体的获取办法可以采用如下方式：.gbff格式中，序列数据以“ORIGIN”开头，并且以“//”符号结束，使用Perl语句截取其中部分，并且还需要去除掉代表序列数量的数值。其次由于DNA是双螺旋结构而序列数据是单向的，所以数据处理第二步要对序列数据进行反向互补链的计算，并且与之前的正向序列数据合并成一个文件。

二、不同k值长度下的未出现k-mer子序列(LAUPs)的计算

获取了数据预处理之后的序列数据后，计算某个k值范围下的物种未出现k-mer子序列的流程如图2所示。首先对同一个物种的某个k值长度进行已出现k-mer结果的计算，在具体实施方式中，可以选取例如k-mer计算软件Jellyfish的软件实现基础计算，首先计算出序列的已出现k-mer结果，Jellyfish将结果存在一个大哈希表中。在此过程中，本发明还提供了一种对原始的Jellyfish方法进行改进后的方法：1.因为在计算未出现的k-mer的过程中，只需要考虑是否出现，而不需要考虑它出现了多少次，所以本发明将Jellyfish软件中的关键参数counting field(出现次数计数，单位为bit)设置为1，并且修改代码，无论排列出现了多少次，都只记为出现一次。2.由于原始的Jellyfish一次只能设定一个k值，本发明对底层算法进行了修改，通过底层算法代码的上述修正设置方式，使得方法支持对同一个数据依次进行多个k值的计算，以此减少数据的载入时间，加快计算速度。

在计算出了物种k-mer的结果之后，依据k的长度，枚举出所有可能的基因序列排列，显然，这些可能总共有4^k个。然后将每一个排列都依次去计算好的已出现k-mer结果中查找，如果已出现k-mer结果中没有找到，将这个序列记为未出现的k-mer排列。

三、相同k值长度下共同未出现k-mer子序列(LAUPs)的计算

由于基因组数据存在一定的测序错误，这可能导致单个物种的计算结果有一定的误差。为了一定程度上消除误差，在计算出了多个物种未出现k-mer排列的结果基础上，对多个物种的未出现k-mer排列进行并集运算，如公式(1)所示，其中CommonLAUPs代表多个物种共同拥有的未出现k-mer子序列，i代表物种的序号，k代表k-mer长度，代表物种i在k值长度下的未出现k-mer排列。N指物种的总数。

四、未出现k-mer子序列(LAUPs)的分析方法

未出现k-mer子序列(LAUPs)的分析方法流程如图3所示，通过计算多个物种的未出现k-mer子序列，可以分析首次出现LAUPs的k值大小，得到一个首次出现的k值阈值范围，以此阈值作为参数，进一步缩小计算的复杂度。在一个具体的实施方式中，可以计算多个具有代表性的物种，并且得到首次出现LAUPs的k的范围为[9,14]。在之后的计算中，可以将计算步骤中的k值上限(如图2中参数N)设置为14，进一步优化计算方法。

其次，在分析基因组数据时，碱基的含量分析是常用的分析方法，计算未出现子序列的GC含量和嘌呤(AG)含量，计算公式如公式2所示。其中，k代表k值的长度，num()函数代表在k长的序列中某种碱基的计数。

在一个具体的实施方式中，未出现子序列的GC含量在60％左右，而AG含量在50％左右。

其次，由于比例不足以说明差异显著性，所以可以使用数学统计分析中的，显著性差异方法，比较未出现子序列与出现的子序列的GC含量与AT含量之间、嘌呤(AG)含量和嘧啶(CT)之间的差异是否显著。在一个具体的实施方式中，比较流程图如图4所示，首先，对含量数据进行正态分布的检验，如果数据符合正态分布，进行方差齐性检验，然后进行T检验，得到一个P值。如果数据不符合正态分布，进行Wilcoxon秩和检验，得到一个P值，最后对P值进行分析，在具体的实施方式中，可以对P值取阈值0.05，P值小于0.05时，认为含量之间是具有显著性差异的。这里需要说明的是，该p值的阈值设定，可以依据所针对的分析对象的差异或者精度需求等进行适应性调整，这些常规的适应性调整均应当视为落入本发明的保护范围之中。

最后，由于计算出的未出现的子序列往往有非常多条，为了快速获取这些序列中可能存在的频繁模式，可以对这些子序列进行Motif发现分析。Motif是指DNA中长度较短，具有保守功能的序列片段，从DNA序列中发现Motif的过程称为Motif发现。在一个具体的实施方式中，可以使用生物信息领域常用的Motif发现工具：MEME，进行Motif发现的基本计算，首先，需要将未出现子序列放在一个文件中，并且满足MEME分析的格式，然后使用MEME在线服务进行Motif发现，发现的结果，是一些Motif Logo图。

实施例2：

在又一个具体的实施例中，本发明还提供了一种多物种未出现k-mer子序列计算和特征分析的系统，所述系统包括：

数据预处理模块，用于获取原始的物种全基因组数据，提取序列数据，并获得反向互补链数据，所述序列数据与反向互补链数据构成预处理数据；

同物种计算模块，用于基于所述预处理数据，进行同一物种不同k值长度下的未出现k-mer子序列计算；

多物种计算模块，用于对多个物种的全基因组数据，进行多物种之间，相同k值长度下共同未出现k-mer子序列的计算；

特性分析模块，用于基于所述步骤2中未出现k-mer子序列以及步骤3中共同未出现的k-mer子序列的数据结果，进行特性分析。

优选的，所述特性分析模块进一步包括：

差异显著性计算单元，用于比较未出现k-mer子序列与出现的k-mer子序列的GC含量与AT含量之间、嘌呤含量和嘧啶之间的差异显著性；

频繁模式计算单元，用于通过Motif发现方式，获取未出现k-mer子序列中可能存在的频繁模式。

需要指出的是，该实施例中的上述系统，其在进行运算时所采用的方法，均可以适用如实施例1中的各种方式。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access Memory，RAM)等。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。