生物序列聚类和全长转录组中转录本亚型识别方法和装置

摘要

本发明公开了一种生物序列聚类和全长转录组测序数据中转录本亚型识别的方法和装置，其中生物序列聚类方法，包括：提供待聚类的生物序列；将所述生物序列比对到参考序列，得到结果文件；按照设定的聚类标准对所述结果文件进行聚类。本发明的生物序列聚类方法，在进行序列相似度计算时加入参考序列信息，能够有效降低计算复杂度，具有快速、准确的特点。

说明书

技术领域

本发明涉及生物序列分析技术领域，尤其涉及一种基于参考序列的生物序列聚类方法和装置以及全长转录组测序数据中转录本亚型识别方法和装置。

背景技术

在生物序列的相似性分析中，相似度的计算可以分为基于序列比对的序列相似度计算方法和基于序列特征分析比较的相似度计算方法，基于序列比对的序列相似度计算方法中采用最多的是序列之间两两比对的策略，基于序列特征分析比较的相似度计算方法主要包括序列的字统计分析、序列之间的编辑距离计算、序列的理化性质统计分析等。

序列相似性分析应用于研究生物的系统进化、同源序列的寻找、大规模生物序列测序的聚类等。例如，在转录组测序的数据分析中，转录本重构是一大难题，相比于二代测序短读取(reads，测序仪输出的核酸序列)拥有大数据量但是很难重构得到全长的信息量，PacBio(Pacific Biosciences公司的简称)测序凭借其超长的读取优势，其Iso-Seq(PacBio的RNA测序商标名)RNA测序几乎可以直测全长转录本，从而能够得到转录本完整的信息量。但是对于PacBio测序，如何从纷繁的测序数据中整理出清晰、完整的信息仍然是很大的挑战，虽然其官方宣称可以直测全长、无需组装进行转录本重构，但是测序数据中多条读取往往拷贝自同一条转录本亚型(isoforms)，从原始下机的全长数据中对转录本亚型进行识别归类仍需要很大的计算量。

目前对于PacBio Iso-Seq测序数据进行转录本重构以及转录本亚型识别的分析只有官方的分析套件(SMRT analysis)，其使用ICE(Isoform level clustering)为核心的算法进行读取聚类(reads clustering)，即转录本亚型识别，ICE算法的描述如下：假如测序得到N(以PacBio RS II测序仪为例，一个测序芯片的N一般为大于一万条序列)条读取，(1)将N条读取用比对软件(如blasr)进行两两比对；(2)然后解析比对结果，标记转录本亚型击中(isoforms hit)关系，建立N*N的转录本亚型击中关系表，一般而言，如果两条序列在比对中没有大的空缺即认为这两条序列来源于同一条转录本亚型；(3)解析转录本亚型击中关系表，形成初步的分组，每个组代表来源于同一条转录本亚型的所有读取拷贝；(4)将各组内的所有读取用pbdagcon(用来对PacBio测序数据构建一致性序列的软件)构建成一条初步的转录本亚型；(5)重新对读取分配，将读取用blasr比对回去上一步中构建的初步的转录本亚型，将读取重新分配给最相似的组，更新步骤(3)中形成的分组，重复步骤(4)；(6)重新对组进行分配，将步骤(5)中得到的转录本亚型互相比对，将相似的组合并到一起，更新步骤(5)中形成的分组，重复步骤(4)；(7)重复步骤(5)、(6)，迭代一定次数后停止。

由此可见，官方的分析套件对于Iso-Seq RNA测序的数据分析是不区分有参考基因组和无参考基因组的，统一用无参基因组的方法进行数据处理，其数据处理过程复杂，耗时耗内存，对于原始数据中包含较高错误率的长读取进行转录本亚型识别需要进行多次比对、聚类、去冗余，计算复杂度很高，随着数据量的增大，计算复杂度的增长也随之增大。

发明内容

本发明提供一种基于参考序列的生物序列聚类方法和装置以及全长转录组测序数据中转录本亚型识别方法和装置，该方法在进行序列相似度计算时加入参考序列信息，能够有效降低计算复杂度，具有快速、准确的特点。

根据本发明的第一方面，本发明提供一种基于参考序列的生物序列聚类方法，包括：提供待聚类的生物序列；将上述生物序列比对到参考序列，得到结果文件；按照设定的聚类标准对上述结果文件进行聚类。

作为进一步改进的方案，上述方法还包括：将上述结果文件分成子文件；相应地，按照设定的聚类标准对上述子文件进行聚类。

作为进一步改进的方案，上述生物序列选自DNA序列、RNA序列或氨基酸序列，上述参考序列选自参考DNA序列、参考RNA序列或参考氨基酸序列。

根据本发明的第二方面，本发明提供一种基于参考序列的生物序列聚类装置，包括：输入单元，用于提供待聚类的生物序列；比对单元，用于将上述生物序列比对到参考序列，得到结果文件；聚类单元，用于按照设定的聚类标准对上述结果文件进行聚类；输出单元，用于输出聚类结果。

根据本发明的第三方面，本发明提供一种基于参考序列的全长转录组测序数据中转录本亚型识别的方法，包括：提供全长转录本序列；将上述全长转录本序列比对到参考序列，得到结果文件；对上述结果文件的序列进行比较外显子结构以及比较每个外显子的5’端和3’端坐标，并且按照设定的标准定义转录本亚型击中；根据上述转录本亚型击中对上述结果文件的序列标记分组；将各组内的所有序列构建成一条转录本亚型的一致性序列。

作为进一步改进的方案，上述方法还包括：将上述结果文件按照设定的标准分成子文件；相应地，对上述子文件的序列进行比较外显子结构以及比较每个外显子的5’端和3’端坐标，并且按照设定的标准定义转录本亚型击中；根据上述转录本亚型击中对上述子文件的序列标记分组；将各组内的所有序列构建成一条转录本亚型。

作为进一步改进的方案，上述外显子结构包括外显子个数以及排列组合方式；上述设定的标准包括：如果两条序列的外显子结构一致，除起始和终止外显子之外的各个外显子容许3bp的错位，则定义为转录本亚型击中；或者，如果两条序列的外显子结构一致，并且对转录本的5’端容错第一预定长度，3’端容错第二预定长度，则定义为转录本亚型击中。

作为进一步改进的方案，上述根据上述转录本亚型击中对上述结果文件的序列标记分组，具体通过基于表搜寻算法的处理来实现。

根据本发明的第四方面，本发明提供一种基于参考序列的全长转录组测序数据中转录本亚型识别的装置，包括：输入单元，用于提供全长转录本序列；比对单元，用于将上述全长转录本序列比对到参考序列，得到结果文件；定义单元，用于对上述结果文件的序列进行比较外显子结构以及比较每个外显子的5’端和3’端坐标，并且按照设定的标准定义转录本亚型击中；分组单元，用于根据上述转录本亚型击中对上述结果文件的序列标记分组；构建单元，用于将各组内的所有序列构建成一条转录本亚型的一致性序列；输出单元，用于输出转录本亚型结果。

本发明的基于参考序列的生物序列聚类方法，在进行序列相似度计算时加入参考序列信息，能够有效降低计算复杂度，具有快速、准确的特点。具体到本发明的基于参考序列的全长转录组测序数据中转录本亚型识别的方法，在有参考序列的情况下，对转录组测序数据分析充分利用参考序列信息，从而能够得到更准确的转录本亚型识别结果。

附图说明

图1为本发明一个实施方案的基于参考序列的生物序列聚类方法的流程图；

图2为本发明一个实施方案的基于参考序列的生物序列聚类装置的结构框图；

图3为本发明一个实施方案的基于参考序列的全长转录组测序数据中转录本亚型识别方法的流程图；

图4为本发明一个实施例中的pass mark算法的示意图；

图5为本发明一个实施方案的基于参考序列的全长转录组测序数据中转录本亚型识别装置的结构框图；

图6为本发明一个实施例中本发明的方法与ICE算法的运行时间比较图。

具体实施方式

下面通过具体实施方式结合附图对本发明作进一步详细说明。

本发明提出一种新的生物序列聚类方法，在进行序列相似度计算时加入参考序列信息，能够有效降低计算复杂度，具有快速、准确的特点。本发明的生物序列聚类方法适用于DNA、RNA和蛋白序列等生物序列的聚类。将需要聚类的生物序列比对到参考序列，比较比对到参考序列的位置从而对序列进行聚类，聚类后的应用有很多，可以去冗余、进行基因表达定量和蛋白表达定量等。

图1示出了本发明一个实施方案的基于参考序列的生物序列聚类方法，包括：

S101：提供待聚类的生物序列。

本发明中待聚类的生物序列可以是DNA、RNA和蛋白序列等，可以来源于测序所得的序列。在本发明的一个实施例中，待聚类的生物序列是全长转录本序列，尤其是来源于PacBio Iso-Seq下机测序数据中的序列，PacBio Iso-Seq下机测序数据可以使用SMRTanalysis分析套件进行去接头、构建插入片段的一致性序列，并从中筛选出三代全长转录本序列。

S102：将上述生物序列比对到参考序列，得到结果文件。

参考序列可以是任何可以用于与待聚类的生物序列进行比对的序列，针对待聚类的生物序列的不同，参考序列也相应的不同，例如，在待聚类的生物序列是DNA、RNA的情况下，参考序列可以是参考基因组序列；在待聚类的生物序列是蛋白序列的情况下，参考序列可以是蛋白家族的序列。在本发明的一个实施例中，待聚类的生物序列是全长转录本序列，参考序列是参考基因组。作为参考基因组，可以是公共的，也可以是由其它途径获得的，只要能够对应上的即可。例如，如果对于人类，可以是NCBI数据库中版本37.3(hg19；NCBIBuild37.3)的人类基因组参考序列。

将上述生物序列比对到参考序列，可以通过目前常用的各种比对方法实现。例如，可以使用GMAP软件将全长转录本序列比对到参考基因组上。为了提升运算速度，在GMAP比对阶段可以对全长转录本序列进行数据分割，就是将整份数据分成N份然后并行处理，增加数据并行程度从而进一步提升运算速度。

结果文件即生物序列比对到参考序列后得到的比对上的文件，其中包含了生物序列在参考序列上的位置等信息，因此可以根据一定的标准对结果文件进行聚类。结果文件的类型可以是各种适用于DNA、RNA和蛋白序列的文件类型，例如常用的gff文件，然而在本发明中，gff文件不是唯一的文件类型，只要能够包含足够信息的文件类型都可以，一般为TAB分隔的文本文件。

S103：按照设定的聚类标准对上述结果文件进行聚类。

聚类标准例如可以是参考序列上的区段范围，例如将规定区段范围内的生物序列定义为一类序列。当然，聚类标准也可以是其它标准，例如对于RNA和蛋白质序列而言，聚类标准可以是基因家族或蛋白家族，例如将属于一个基因家族或蛋白家族的RNA或蛋白质序列定义为一类序列。

在本发明的一个实施例中，在待聚类的生物序列是全长转录本序列，参考序列是参考基因组的情况下，聚类标准可以是外显子结构(例如，外显子个数以及排列组合方式)以及每个外显子的5’端(核酸序列的起始端)和3’端(核酸序列的终止端)的坐标，即根据上述外显子结构和坐标而将全长转录本序列进行聚类。

本发明的方法可以直接对结果文件进行聚类，然而，在某些情况下，结果文件可能数据量非常巨大，直接对结果文件进行聚类，可能产生大量的冗余计算量，从而导致运算变慢。因此，作为进一步改进的方案，可以将结果文件分成子文件，也就是将大的计算问题拆分成小块，减少冗余计算量；然后，按照设定的聚类标准对子文件进行聚类。

对应于上述基于参考序列的生物序列聚类方法，本发明还提供一种基于参考序列的生物序列聚类装置，如图2所示，该装置包括：输入单元201，用于提供待聚类的生物序列；比对单元202，用于将生物序列比对到参考序列，得到结果文件；聚类单元203，用于按照设定的聚类标准对结果文件进行聚类；输出单元204，用于输出聚类结果。

作为本发明的基于参考序列的生物序列聚类方法的一个具体应用，本发明还提供一种基于参考序列的全长转录组测序数据中转录本亚型识别的方法，在有参考序列(例如，参考基因组)的情况下，例如对PacBio转录组数据分析充分利用参考基因组信息，从而能够得到更准确的转录本亚型识别结果，基于PacBio全长转录本在转录本亚型识别上的计算量大、软件处理流程复杂的不足，本发明所提供的转录本亚型识别的方法能够有效降低计算复杂度，快速，准确。

如图3所示，一种基于参考序列的全长转录组测序数据中转录本亚型识别的方法包括：

S301：提供全长转录本序列。

全长转录本序列，尤其是来源于PacBio Iso-Seq下机测序数据中的序列，PacBioIso-Seq下机测序数据可以使用SMRT analysis分析套件进行去接头、构建插入片段的一致性序列，并从中筛选出三代全长转录本序列。

需要说明的是，本发明所述的全长转录本为绝对定义，并非PacBio所特有，也即以其它技术获得全长转录本，然后进行基于参考序列比对进行外显子结构解析的转录本亚型识别方法都应为本发明所包含的范畴。

S302：将上述全长转录本序列比对到参考序列，得到结果文件。

在本发明的一个实施例中，待聚类的生物序列是全长转录本序列，参考序列是参考基因组。作为参考基因组，可以是公共的，也可以是由其它途径获得的，只要能够对应上的即可。例如，如果对于人类，可以是NCBI数据库中版本37.3(hg19；NCBI Build 37.3)的人类基因组参考序列。

作为一种变形，可以使用非基因组的其它参考序列进行比对然后解析外显子结构，例如用全长转录本测序数据构建出每个基因家族的参考序列，然后再将全长转录本比对回此参考序列。因此，本发明所述的参考序列是包含参考基因组在内的其它任何能够用来作为全长转录本采用比对的形式进行外显子结构解析进而进行转录本亚型识别的参考序列。

将上述生物序列比对到参考序列，可以通过目前常用的各种比对方法实现。例如，可以使用GMAP软件将全长转录本序列比对到参考基因组上。为了提升运算速度，在GMAP比对阶段可以对全长转录本序列进行数据分割，就是将整份数据分成N份然后并行处理，增加数据并行程度从而进一步提升运算速度。

结果文件即全长转录本序列比对到参考序列后得到的比对上的文件，其中包含了全长转录本序列在参考序列上的位置等信息，因此可以根据一定的标准对结果文件进行聚类。结果文件的类型可以是各种可用的文件类型，例如常用的gff文件，然而在本发明中，gff文件不是唯一的文件类型，只要能够包含足够信息的文件类型都可以，一般为TAB分隔的文本文件。

本发明的方法可以直接对结果文件进行聚类，然而，在某些情况下，结果文件可能数据量非常巨大，直接对结果文件进行聚类，可能产生大量的冗余计算量，从而导致运算变慢。因此，作为进一步改进的方案，可以将结果文件按照设定的标准(例如染色体、正负链等)分成子文件，也就是将大的计算问题拆分成小块，减少冗余计算量；然后，按照设定的聚类标准对子文件进行聚类。例如，在本发明的一个实施例中，将gff结果文件根据染色体和正负链的关系分割成2X(X为染色体个数)个gff子文件，相应地，后续定义转录本亚型击中、分组和构建转录本亚型的步骤，也是针对子文件进行的，以下仅以对结果文件的处理为例加以说明，在分成子文件的情况下，对子文件的处理是相同的。

S303：对上述结果文件的序列进行比较外显子结构以及比较每个外显子的5’端和3’端坐标，并且按照设定的标准定义转录本亚型击中。

外显子结构包括外显子个数以及排列组合方式等，外显子的5’端即核酸序列的起始端，外显子的3’端即核酸序列的终止端。设定的标准，例如可以是：如果两条序列的外显子结构一致，除起始和终止外显子之外的各个外显子容许3bp的错位，则定义为转录本亚型击中；或者，如果两条序列的外显子结构一致，并且对转录本的5’端容错第一预定长度，3’端容错第二预定长度，则定义为转录本亚型击中。

上述标准考虑到实际情况，引入了容错机制。具体地，考虑到外显子容易存在3bp的错位的情况，以及转录本的5’端和3’端经常会发生降解等情况。第一预定长度和第二预定长度可以根据具体的实验情况来确定，一般而言，在真实情况下，在Iso-Seq文库构建过程中，转录本的5’端更容易降解，3’端由于有polyA(连续的A碱基组成的序列)序列保护，转录本的3’端在基因组上的坐标比5’端在基因组上的坐标更加准确，因此本发明可以对转录本的5’端容错50bp左右，3’端容错30bp左右，达到条件则定义为转录本亚型击中。

现有的ICE算法不利用参考序列信息构建转录本亚型，最终比回参考序列会造成一定的误差，相比而言，本发明的方法，一方面加入参考序列的信息，对转录本亚型边界能定义得更加准确，另一方面采取了符合实际的容错机制，更加准确。

此外，对转录本起始端坐标容错区间的优化，可以根据测序数据的坐标的分布趋势定义更加合理的容错区间。

S304：根据上述转录本亚型击中对上述结果文件的序列标记分组。

分组可以采用一种基于表搜寻算法(pass mark算法)的处理来进行。算法描述如下：(1)对每个结果文件或gff子文件包含的N条读取构建N*N的循环进行比较；(2)遍历循环的第一行，根据转录本亚型击中标记分组，标记已有属组的读取；(3)跳到还没有分配到属组的读取开头的行，根据转录本亚型击中标记分组，一行中已有属组的读取将不进行计算转录本亚型击中；(4)重复步骤(3)，直到所有读取都标记上了分组。

图4示出了一种示例性的pass mark算法的图示，这里假设有ABCDE五条序列需要分组，表中每个格子用数字1表示两条序列为转录本亚型击中关系，数字0表示非转录本亚型击中关系，从图上可以明显的看出pass mark算法只需要计算整个表25个格子中的5个(灰色填充的格子)就可以将五条序列进行分组。由此可见，比对次数大大减少，构建转录本亚型击中的步骤只需要比对N次，而ICE需要比对N*N次，对转录本亚型击中的关系表能够根据染色体和正负链进行拆分，从而缩减了表的大小，减少无关计算，在寻找分组的过程中采用了pass mark算法，进一步压缩了分析时间，数据处理步骤更加简便，无需反复迭代过程。

S305：将各组内的所有序列构建成一条转录本亚型的一致性序列。

可以使用pbdagcon软件(用来对PacBio测序数据构建一致性序列的软件)将各组内的所有读取构建成一条转录本亚型。

本发明提出方法与ICE算法的运行时间比较

对应于本发明的基于参考序列的全长转录组测序数据中转录本亚型识别的方法，本发明还提出一种基于参考序列的全长转录组测序数据中转录本亚型识别的装置，如图5所示，该装置包括：输入单元501，用于提供全长转录本序列；比对单元502，用于将上述全长转录本序列比对到参考序列，得到结果文件；定义单元503，用于对上述结果文件的序列进行比较外显子结构以及比较每个外显子的5’端和3’端坐标，并且按照设定的标准定义转录本亚型击中；分组单元504，用于根据上述转录本亚型击中对上述结果文件的序列标记分组；构建单元505，用于将各组内的所有序列构建成一条转录本亚型；输出单元506，用于输出转录本亚型结果。

本领域普通技术人员可以理解，上述实施方式中各种方法的全部或部分步骤可以通过程序来指令相关硬件完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

以下通过实施例详细说明本发明的具体实现和技术效果，应当理解，实施例仅是示例性的，不能理解为对本发明保护范围的限制。

实施例

本实施例选取PacBio公开的MCF7数据集2015批次的测序数据为测试数据进行算法的性能上的比较，数据来源于：http://datasets.pacb.com.s3.amazonaws.com/2015/IsoSeqHumanMCF7Transcriptome/list.html。

本发明的方法和用于对比的ICE算法的计算机软、硬件运行环境均是：软件为CentOS系统、硬件为双核16线程、64G内存。

本发明的方法的具体实现步骤如下：

1、使用SMRT analysis分析套件对PacBio Iso-Seq下机测序数据进行去接头、构建插入片段的一致性序列，并从中筛选出三代全长转录本序列；

2、将全长转录本序列比对回参考基因组，得到比对结果的gff文件；

3、将gff文件根据染色体和正负链的关系分割成2X(X为染色体个数)个gff子文件；

4、对每个gff子文件的序列进行比较外显子结构(外显子个数以及排列组合方式)以及比较每个外显子的5’端(核酸序列的起始端)和3’端(核酸序列的终止端)坐标，如果两条序列的外显子结构以及起始和终止坐标一致(除开起始和终止外显子，容许3bp的错位)则定义为转录本亚型击中；

5、对转录本的5’端容错50bp，3’端容错30bp，达到条件则定义为转录本亚型击中；

6、对每个gff子文件进行pass mark算法(一种基于表搜寻算法)处理；

7、将各组内的所有reads用pbdacon构建成一条转录本亚型。

表1示出了本发明的方法与ICE算法的运行时间比较结果。

表1本发明的方法与ICE算法的运行时间比较

表1的结果显示：不同芯片数量的测序数据对两种算法的运行时间进行比较，可以明显看到本发明的方法在分析时间上具有明显的优势。

图6示出了本发明的方法与ICE算法的运行时间比较图。从图可以看出，ICE算法在数据量增大的时候在分析时间上具有更加陡峭的增长趋势，而本发明的方法则对数据量的增加表现出更平缓的增长趋势，说明本发明的方法能够很有效的降低分析复杂性。

表2示出了本发明的方法与ICE算法的转录本亚型识别准确度的比较结果。

表2本发明的方法与ICE算法的转录本亚型识别准确度的比较

注：“校正”为选用测序数据中的非全长转录本对识别出来的转录本亚型进行测序错误的校正，校正选用官方提供的quiver程序进行；“替换”、“删除”、“插入”为基于序列比对的三种最基本的单碱基突变，这里将校正前后的转录本亚型比对回参考基因组，然后统计其三种突变比例，比例越低说明转录本亚型越准确；“高质量转录本亚型”为识别出来的转录本亚型质量值达到99％以上的准确率，也即100个碱基中有小于或者等于1个碱基错误，在PacBio Iso-Seq分析中最终进行下游结构与功能分析采用的就是高质量转录本亚型部分，所以高质量转录本亚型数量越多说明方法越有优势；“一致性”为两种方法识别出来的转录本亚型进行互相比对，统计具有转录本亚型击中关系的转录本亚型占识别到的转录本亚型数的比例，一致性用来说明两种方法的基本框架是否能一致。

表2的结果显示，在选取的测试数据上，本发明的方法在校正前后的三种单碱基突变比例都比ICE算法更低，说明本发明的方法所识别的转录本亚型更加准确；可以看到采用容错机制后识别到的转录本亚型数上，本发明的方法比ICE要少；可以看到本发明的方法在识别到的高质量转录本亚型数量更多，这将能够提供更加好的数据利用率，保留了更多的可用信息；可以看到一致性方面，两种方法识别的转录本亚型有转录本亚型击中关系的序列都大于99％，说明两种方法在容易识别的部分是一致的，区别在于本发明的方法有更好的容错机制，更加合理，从而也更准确。

以上内容是结合具体的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。