一种基于序列模式挖掘算法的系统发生树构建方法与流程

本发明涉及一种系统发生树构建方法，具体涉及一种基于序列模式挖掘算法的系统发生树构建方法。

背景技术：

自上世纪八十年代开始，伴随计算机技术，测序技术和分子生物学的不断发展，各国研究人员实施完成了多项基因组计划，通过基因组测序、蛋白质测序和结构解析等方法，积累了大量分子生物学数据。然而面对这些以指数速率增长的数据，研究者却只获得了少量有价值的信息，大量具有生物学意义的潜在信息淹没在各种各样的数据库中，人们分析处理数据的能力已经远落后于数据产生的能力。这就迫切需要计算能力更强大，硬件价格更优惠的计算设备以及数据分析研究者开发出复杂度更低、结果更优的算法来帮助人们分析处理各种各样的测序数据。面对种类繁多的测序数据，快速，高效地挖掘具有生物学意义的序列信息，已经是一个科学研究的热点问题，其中通过序列比对的方法，确定dna、rna、蛋白质序列之间可能的结构、功能、进化上的相似性是生物信息学中一个重要的基础研究。通过将未知序列和已知序列进行相似性比较，可以推测未知序列所具有的功能和特点，因为相似的序列很可能具有相似的结构，相似的结构也很可能具有相似的功能。以进化学说为理论基础，借助序列比对结果，如果比较的两条序列之间相似性足够高，那么就可推断这两条序列从同一个祖先进化而来，在进化分析中意义重大。

分子生物学中构建系统发生树一个重要且关键的步骤就是对序列进行相似性比较。衡量序列之间的相似性关系，目前的技术策略主要可分为基于多序列比对的方法和非序列比对的方法。

1，基于多序列比对的方法。

比较两条或是多条序列之间是否相似，传统方法一般采用基于动态规划算法的方法。简单来说，多序列比对的过程就是通过在序列中插入空格的方式，来让序列中更多相似的区域对齐的过程。这在序列相似度较高，序列集较小的时候，通常可以得的一个可靠的比对结果。但随着分子生物学的快速发展，人们测序积累了大量不同物种不同器官组织的序列数据，一次比对任务，参与比对的序列数目成百上千，且序列间的相似性可能很低，这就导致传统的比对方法，无法得到一个可靠的比对结果，或是直接无法比对。代表性的多序列比对工具有clustalw家族软件和muscle等。然而随着人们对复杂进化途径的理解和关于生物序列性质的熟悉，仅仅采用基于多序列比对的序列比较方法的缺点逐渐显露。主要有以下几点：

(1)基于比对的序列比较方法假设同源序列包含一系列线性排列和或多或少保守序列的延伸，然而这些假设通常与现实情况相违背。

(2)进化过程中，随着序列不断积累单碱基突变和一些小片段的插入和删除，序列之间的相似度越来越低，当序列相似性低于一定临界点时，基于比对方法的准确性将迅速降低，甚至无法比对。

(3)由于大多数多序列比对方法基于动态规划算法，时间复杂度高，资源占用多，尤其当有大规模低相似度的序列需要比对时，将是一个耗时耗力的过程。同时多序列比对计算是一个np-hard问题。

(4)最终的比对结果会受到多个因素的影响，比对指导树的构建，替换计分矩阵的选择，空位罚分机制的设定等一些先验假设。

(5)两条序列之间的距离过度依赖于那些已经对齐的区域，忽略了一些本身含有一定生物学意义的序列片段。

2，非比对的方法。

非比对方法的出现，一定程度上克服了多序列比对的一些缺点，这类方法的应用已经从全基因组种系分析发展到蛋白质家族分类、序列重排检测等。这些强有力的方法在分析和处理下一代测序数据方面尤为有用。序列比对算法的目的就是尽可能多的定位那些在生物序列中可能具有结构、功能、进化相关的相似性区域，非比对序列比较是一种在量化序列之间相似性的计算过程中，任何一步都不需要将序列对齐的方法。这种方法不依赖于动态规划算法，计算复杂度底，不依赖于序列进化过程中的种种假设，适合全基因组比较，同时在序列保守区域较少和存在序列重排等情况下仍可工作。目前常用的非比对的方法除了部分基于信息论和图的方法，主要包括两大类：

第一类基于kmer的方法。比如featurefrequencyprofile(ffp)方法和组合向量表示cvtree方法，这些方法使用一个提前设定好大小的滑动窗口在序列集中提取特征。然后将每一条序列转换为特征频率向量的形式，利用这些向量集就可以来计算序列之间的距离。不同与这两种方法使用精确匹配模式的方法来提取特征，spacedword的方法通过在模式中加入模糊匹配位的方式对模式进行扩展，然后计算序列之间的距离。但是采用这种策略的方法，有一个不可避免的缺点。这些方法在提取特征的时候提前设定了模式的大小，这在序列相似度较低，有大规模序列需要比较时往往会提取一些干扰模式降低后续计算速度，占用大量资源。

第二类基于substring的方法。比如averagecommonsubstring(acs)方法和underlyingapproach(ua)方法。相比于基于kmer的方法，基于substring的方法提取模式的长度是可变的，用户不用自己提前设定模式大小。但是他们在计算每两条序列之间的距离时，仅仅利用当前两条序列之间的共有模式来计算序列之间的相似度，忽略了全局序列之间潜在的关系，忽略了每个序列子集之前的区别与联系。

技术实现要素：

本发明的目的在于提供一种快速准确、能够对大规模、低相似度的基因组数据进行系统发生树构建的基于序列模式挖掘算法的系统发生树构建方法。

为实现上述目的，本发明是通过以下技术方案来实现：

一种基于序列模式挖掘算法的系统发生树构建方法，包括如下步骤：

步骤1：挖掘隐藏在序列集中能够用于衡量序列相似性的频繁序列模式，得到初始模式集；

步骤2：过滤初始模式集中的非闭合频繁模式，得到能够代表序列集的优化模式集；然后构造模式向量集，再计算数值向量之间的距离，进而构建距离矩阵；根据距离矩阵生成最终需要的系统发生树。

本发明进一步的改进在于，步骤1的具体步骤如下：

第一，根据输入序列集大小、序列类别以及序列长度信息，设定模式挖掘过程中支持度参数、精确匹配与模糊匹配位参数，指定输入文件名及路径、输出特征集文件名及路径以及输出距离矩阵文件名及路径；

第二，根据设定的支持度参数、精确匹配与模糊匹配位参数，通过序列模式挖掘算法对序列集进行挖掘，得到初始模式集。

本发明进一步的改进在于，通过序列模式挖掘算法对序列集进行挖掘的具体过程为：初始模式以一个空模式和其相对应的空投影数据库开始模式增长过程，在模式增长的过程中，陆续保留那些满足支持度参数、精确匹配与模糊匹配位参数的频繁模式；相应的投影数据库迭代更新，直到投影数据库无法更新以满足支持度参数，模式停止增长，得到初始模式集。

本发明进一步的改进在于，步骤2中，构造模式向量集时：将初始模式集中每条序列转换为二进制的模式向量的形式，或是利用模式自身所含信息量的多少来将序列转换为赋有权重信息的模式向量的形式。

本发明进一步的改进在于，步骤2中，构造模式向量集时，若构造赋有权重信息的模式向量，采用利用模式的长度信息的方式、去除模式之间的重叠部分的方式、对模式进行合并及扩展操作的方式、组合具有先后顺序的模式为大模式的方式中的一种或几种构造赋有权重信息的模式向量。

本发明进一步的改进在于，采用利用模式的长度信息的方式构造赋有权重信息的模式向量时，权重通过以下公式计算得到；

其中，w为权重，patlength表示当前模式的长度，dbsize和patsup分别表示序列数据库的大小和当前模式的支持度大小。

与现有技术相比，本发明具有的有益效果：

本发明以快速准确构建系统发生树为出发点，采取序列模式挖掘算法的方法来提取序列集中频繁出现的特异性模式。然后利用过滤后的模式集将序列集转换为模式向量集的形式计算得到距离矩阵，最后利用距离矩阵构建系统发生树。针对大规模，低相似度序列集的特点，基于序列模式挖掘算法的方案，利用模式增长的策略可以挖掘序列集中更具代表性的模式，避免了提取对衡量序列相似性无用的冗余模式，优化了全局范围内序列之间相似性的衡量。其具有的优点如下：

第一：本发明是一种非比对的系统发生树构建方法，相比较于传统的基于多序列比对的方法，本发明，可以快速构建可靠的系统发生树。

第二：由于序列不断积累单碱基突变和小片段插入和删除等变异，序列之间的相似性变得越来越低。对于低相似度的序列集，采用模式挖掘的方法可以提取个别序列子集中包含的序列片段，更好的对序列集进行分类。

第三：大规模序列集对基于多序列比对的系统发生树构建方法来说是一个巨大的挑战，往往需要耗费大量的时间来进行序列的比对任务，结果还可能是一个不可靠的比对结果。使用模式增长的方法提取序列中的特异性模式，往往会在接近线性的时间内得到一个满意的模式集。利用模式集可快速的构建一颗可靠的系统发生树。

第四：按需设定好初始模式挖掘算法参数，就可直接通过模式增长的方法得到特征集。不需要像多序列比对的方法，依靠假设，设定空位罚分，选择替换计分矩阵等进行多序列比对过程。

第五：基于多序列比对的方法，过分强调序列中已经对齐的区域，忽略了个别序列中包含的具有生物学意义的序列片段。本发明不在强调对齐区域，只要某序列片段频繁出现在多条序列中，满足设定的最小初始条件，就会被用于后续的进化树构建过程。

第六：区别于基于kmer的方法，使用一个固定大小的滑动窗口来提取序列集中的特征，基于模式挖掘方法的特征提取过程是一个长度非依赖的过程。频繁模式的长度在满足精确匹配位数和模糊匹配位数的要求下，闭合的频繁序列模式都会保存到序列模式集中。尤其是对于低相似度的序列集来说，如果提取特征的时提前设定了模式的大小，有多条序列需要比较时往往会提取一些干扰模式降低计算速度，占用大量资源。

第七：区别于基于substring的方法，特征的提取只在两条序列之间得到，本发明以全部序列为对象进行频繁模式的挖掘。

第八：计算全局范围内序列之间的相似性，是本发明优于substring或是采用lempel-ziv复杂度策略的非比对方法的一个方面。

附图说明

图1为现有的基于多序列比对的系统发生树构建大致流程图。

图2为本发明的基于序列模式挖掘算法的系统发生树构建流程图。

图3为对比不同方法在序列相似度较低时的比较结果。

图4为对比不同方法在不同大小序列集上的比较结果。

具体实施方式

下面结合附图对本发明进行详细描述。

本发明的一种基于序列模式挖掘算法的系统发生树构建方法，包括如下2个步骤。

步骤1：为了准确地计算序列之间的相似性，克服基于多序列比对方法的复杂计算过程与各种推测假设，克服基于kmer或是substring等方法提取过多干扰模式与忽略全局范围内序列之间的相似性关系的不足，开发基于序列模式挖掘算法的方案来对大规模，低相似度的序列集，挖掘频繁出现在多条序列中的特异模式，然后基于这些模式来构建系统发生树是本发明的目的。步骤1所完成的任务的就是挖掘隐藏在序列集中能够用于衡量序列相似性的特异模式。挖掘频繁序列模式的具体步骤如下：

主要包括二个步骤。第一步，根据输入序列集大小(有多少条序列)，序列类别(dna、蛋白质序列)，序列长度信息，设定模式挖掘过程中支持度参数(最少有多少条序列含有当前模式)，模式中精确匹配位(与序列中相应位置的符号信息一一对应的位数)和允许连续出现的模糊匹配位(可以连续匹配序列中任意符号的位数)大小，指定输入文件名及路劲，输出特征集文件名及路径，输出距离矩阵文件名及路径。

第二步，根据按需求设定的支持度参数、精确匹配与模糊匹配位参数，通过挖掘算法对序列集进行挖掘，得到初始模式集。

通过序列模式挖掘算法对序列集进行挖掘的具体过程为：初始模式以一个空模式和其相对应的空投影数据库开始模式增长过程，即模式以一个空根开始增长。在模式增长的过程中，陆续保留那些满足支持度参数、精确匹配与模糊匹配位参数的频繁模式；相应的投影数据库迭代更新，直到投影数据库无法更新以满足支持度参数，模式停止增长，得到初始模式集。

步骤2：满足步骤1中支持度参数、精确匹配与模糊匹配位参数的初始模式集中会有部分模式是非闭合的频繁模式，模式之间也会有重叠存在，所以在得到初始模式集后可以过滤掉非闭合频繁模式，得到更能代表序列集的优化模式集。

利用优化模式集来衡量序列之间的相似性，所采取的方法是将序列集中每条序列转换为二进制的模式向量的形式，或是利用模式自身所含信息量的多少来将序列转换为赋有权重信息的模式向量的形式。利用转换得到的模式向量集，就可以计算数值向量之间的距离，进而构建距离矩阵。根据距离矩阵就可以利用现有的一些方法来生成最终需要的系统发生树。

构造模式向量集时，若构造赋有权重信息的模式向量，采用利用模式的长度信息的方式、去除模式之间的重叠部分的方式、对模式进行合并及扩展操作的方式、组合具有先后顺序的模式为大模式的方式中的一种或几种构造赋有权重信息的模式向量。

本发明所述的基于序列模式挖掘算法的系统发生树构建，如图2所示，包括两大核心模块，分别为序列模式挖掘模块，模式过滤、去重、组合、赋权重等模块。同时包括利用优化后的模式集来计算距离矩阵模块，利用最近邻等方法来构建系统发生树模块。两大核心模块具体详述如下。

1序列模式挖掘模块。

相比于图1所示现有的经典经过多序列比对后利用最大似然法等构建系统发生树的方法，序列模式挖掘方法的模式提取是一个非比对的过程。对比于基于kmer或substring等方法的非比对系统发生树构建方法，提取可变长度的频繁模式与衡量序列间全局相似性，基于序列模式挖掘的方法有着明显的优势。以大规模，低相似度的序列集举例来说，模式挖掘主要分二个步骤：

第一，对于输入序列集，如利用rose(generatingsequencefamilies)软件生成不同相似度，不同数目的fasta格式的序列集。首先将序列名称和相对应的序列存储于序列标签集和序列数据集中，规定挖掘得到的频繁模式包含其支持度信息和位置信息。设定初始挖掘条件，包括支持度大小，精确匹配位和连续出现的最长模糊匹配位大小，指定模式集和距离矩阵输出文件名及地址信息。

第二，模式从一个空模式和其相应的空投影数据库开始增长，以dna为例，首先模式从空模式增长变为一个碱基长度的模式，相应的投影数据库更新为序列集中所有当前碱基在哪一条序列中哪些位置出现的信息，此时模式完成一次增长。模式继续增长，从之前一个碱基长度的模式增长为两个碱基长度的模式，检测投影数据库中每一个位置的下一个位置是否是新增长的碱基。如果是，更新当前的位置信息为下一个碱基出现的位置。如果不是，删除当前位置。投影数据库更新过程中，不断删除那些不包含任何位置信息的序列。如此，模式不断增长，相应的投影数据库不断更新，直到当前的模式不满足设定的初始条件，模式停止增长。在模式增长的过程中，模式新增长的一位可以是一个确切的碱基，也可以是一个匹配任何碱基的模糊匹配位，当模式新增长的一位是模糊匹配位时，相应的投影数据库中所有位置信息都更新为向后加一位。如果向后加一位超出序列长度就将其删除。

输出需要保存的模式之前，检测此模式的支持度和精确匹配位数是否满足设定最小的支持度大小和精确匹配的位数大小，同时检测模式的最后一位是否为模糊匹配位，当模式最后一位是非模糊匹配位，且支持度和精确匹配位数满足设定要求，保存当前模式。

2模式过滤、去重、组合、赋权重等模块。

通过序列模式挖掘模块得到的初始模式集需要在构造模式向量集之前进一步的处理，如过滤非闭合模式，或附加去除模式之间的重叠部分，对模式进行扩展，给模式赋予一定的权重，将具有先后顺序的模式组合为一个更大的模式。有利于提高后续计算速度和构建的系统发生树的质量。

当一个模式的子模式和其具有相同的支持度时，这个模式被认为是非闭合的模式。通过上述模式增长的方式得到的初始模式会包含非闭合的频繁模式，为了减小得到的模式集的大小，加快距离矩阵的计算速度，在模式增长的过程中添加闭模式检测环节。检测的方法是在模式更新的同时，检测模式增长前的模式是否与模式增长后的模式具有相同的支持度，并且模式增长前的模式是模式增长后的模式的子集，如果这两个模式具有相同的支持度，那么删除之前的模式，仅保留模式增长后的模式用于后续模式向量的构建。在去除了非闭合模式后的模式集中，有部分模式之间存在重叠部分，在后续给每个模式按照其长度信息赋予一定权值的时候，如果a模式的一部分是b模式的一部分且b模式已经按其本身长度信息赋予了一定权值，后续给a模式赋权重的时候，为了避免重复计算重叠部分的模式长度，给a模式只赋值其自身与b模式非重叠部分模式长度相关的权值。

按照模式本身组成元素性质的相似性关系，可以将模式进行一定程度的合并与扩展。如模式abc中a的生化性质与d相似，那么模式abc就可能与模式dbc相似，可以将模式abc合并为(a/d)bc。同时生物序列中，某两个序列片段的位置关系经常在序列中含有一定的顺序关系，因此可以将在顺序上存在先后关系的模式链接或扩展为超模式，当序列之间含有这种超模式的时候，认为这两条序列更有可能相似。

共享的序列片段长度越长，序列间越有可能相似。因此不同长度的序列片段应当区别对待，给每个不同长度的序列片段赋予一定的权重，然后根据每个模式的权值来构造计算距离矩阵所用的模式向量。给模式赋的权重公式如下：

其中，w为权重，patlength表示当前模式的长度，dbsize和patsup分别表示序列数据库的大小和当前模式的支持度大小。

利用优化赋权值后的模式集，可以将序列集转换为模式向量集的形式。对于相似度较高的序列集，可将每一条序列转换为关于模式集的二进制模式向量形式，0表示此序列不包含这个模式，1表示此序列含有这个模式，就可得到理想的系统发生树。对于低相似度序列集，将序列转换为模式向量时，若序列包含某个模式，就用模式的权值代替1，序列不包含某个模式，则继续用0表示。

计算数值向量之间的距离，可以采用但不限于jensen-shannon距离、余弦距离、杰卡德距离等。

有了每两条序列之间的距离值，就可以得到关于此序列集的距离矩阵，利用距离矩阵作为输入就可通过最近邻方法来构建一颗系统发生树。

利用rose(generatingsequencefamilies)生成200条和600条长度在350个氨基酸左右，序列之间的相似度小于10％的序列集上做测试。比较不同方法构建的系统发生树与标准参考树之间的robinson-fouldes距离，距离值越小表示与标准参考树越相似。并且测试各种方法在相同硬件条件下的运行速度。

如图3所示，在序列之间的相似度小于10％时，对200条模式序列构建系统发生树。由于序列之间的相似度过低，基于多序列比对的方法无法比对得到可靠的比对结果或是直接无法比对，用于比较的五种方法都是非比对的系统发生树构建方法。分别是featurefrequencyprofile(ffp)方法，组合向量表示cvtree(cv)方法，spacedword(sw)方法，underlyingapproach(ua)方法，averagecommonsubstring(acs)方法和本发明基于序列模式挖掘算法的系统发生树构建方法(spv)方法。五种比对方法的参数使用尽可能得到最好建树结果的参数，基于序列模式挖掘算法的spv方法使用支持度大小为0.03，精确匹配位和模糊匹配位大小分别为2和2。从图中可见，基于序列模式挖掘方法的系统发生树构建方法在序列相似度较低可以得到比其他非比对方法更好的结果。如图4所示，在序列集较大有600条序列需要构建系统发生树时，相比于其他方法，本发明同样可以构建更优的系统发生树。

表1不同方法运行时间比较

如表1所示，第一列表示不同的非比对系统发生树构建方法，其中各个方法后括号中的数值表示建树时使用的参数，若没有则使用推荐的参数值。spv后括号中2-1表示精确匹配位和允许最长模式匹配位数的大小分别为2和1。表中第二列表示各方法构建的系统发生树与标准参数树之间的robinson-fouldes距离，值越小表示与标准参考树越相似。表中第三列表示相同硬件条件下各方法运行时间。从表1中可以发现，在较短的运行时间内基于序列模式挖掘算法的系统发生树构建方法可以得到一个比其他方法更好的结果。