一种基于后缀树算法的基因组测序序列与参考基因组比对的方法

摘要

本发明涉及生物信息技术领域，具体涉及一种基于后缀树算法的基因组测序序列与参考基因组比对的方法。本发明提供的基于后缀树算法的基因组测序序列与参考基因组比对的方法包括构建参考基因组索引以及将基因组测序序列与参考基因组索引进行序列比对的步骤，其中，所述构建参考基因组索引包括如下步骤：(1)构建参考基因组索引的初步后缀树；(2)将所述初步后缀树中含有分叉的节点转换成节点数字，不含有分叉的节点转换成节点矩阵，构建后续用于比对的最终后缀树。本发明提供了一种占用内存相对较小、运行速度较快的、基于后缀树算法进行序列比对的方法，有效降低了读入索引对计算机内存的要求。

说明书

技术领域

本发明涉及生物信息技术领域，具体涉及一种基于后缀树算法的基因组测序序列与参考基因组比对的方法。

背景技术

序列比对是生物信息学的基本组成和重要基础。序列比对的基本思想是：基于生物学中序列决定结构、结构决定功能的普遍规律，将核酸序列和蛋白质一级结构上的序列都看成由基本字符组成的字符串，检测序列之间的相似性，发现生物序列中的功能、结构和进化的信息。

目前主要的序列比对算法分为两种：基于哈希表数据结构比对算法和BWT索引数据结构比对算法。其中，基于哈希表数据结构比对算法是将参考基因组分割成短序列，并将其存入哈希表(一种数据结构)中。在进行序列比对时，通过检索短序列数据集的哈希表，然后找到ref跟reads完全匹配(match)的子序列，即seed位点，而后通过经典的全局比对或者局部比对，将整条reads比对上去。后缀树算法是基于哈希表数据结构比对算法中的一种，具有运行速度快的特点，但比对时需要读入内存的参考基因组构建的索引非常大，使得基于后缀树算法的比对软件通常占用内存过大，对计算机的配置要求很高。

发明内容

本发明的目的是提供一种基于后缀树算法的基因组测序序列与参考基因组比对的方法，该方法将参考基因组构建索引的序列转化成数字，在保留了基于后缀树算法比对软件的运行速度快的优势的同时，极大地降低了传统后缀树算法比对软件占用的内存。

具体地，本发明提供以下技术方案：

首先，本发明提供一种基于后缀树算法的基因组测序序列与参考基因组比对的方法，该方法包括构建参考基因组索引以及将基因组测序序列与参考基因组索引进行序列比对的步骤，

其中，所述构建参考基因组索引包括如下步骤：

(1)构建参考基因组索引的初步后缀树；

(2)将所述初步后缀树中含有分叉的节点转换成节点数字，不含有分叉的节点转换成节点矩阵，构建后续用于比对的最终后缀树。

上述步骤(1)中，所述构建参考基因组索引的后缀树包括：

①将参考基因组序列按照100-200碱基长度、1个碱基的步长进行截断，在每条截断的序列的最后插入$符号，并找到每条截断的序列的所有后缀，将各后缀按照每条后缀第一个字符在原序列中位于的坐标进行编号；

②将所有的后缀按照字母表的顺序进行排序；

③将经过排序的每一个后缀作为节点，按照有序相邻的后缀最长公共前缀将节点合并，构成初步后缀树，每个节点记录字符串的位置。

具体地，以序列AGTCGACAG为例，按照上述方法构建初步后缀树：

①在序列的最后插入“$”符号，并找出序列的所有后缀；

1AGTCGACAG$

2GTCGACAG$

3TCGACAG$

4CGACAG$

5GACAG$

6ACAG$

7CAG$

8AG$

9G$

10$

②将所有的后缀按照字母表的顺序进行排序；

10$

6ACAG$

1AGTCGACAG$

8AG$

7CAG$

4CGACAG$

5GACAG$

2GTCGACAG$

9G$

3TCGACAG$

③构建初步后缀树；

将每一个后缀作为节点，并将独立的节点合并，构成后缀树，每个节点记录字符串的位置(如图1所示)

A G T C G A C A G$

1 2 3 4 5 6 7 8 9 10

上述步骤(2)中，所述最终后缀树的构建包括：

①对于所述初步后缀树中的每条序列，分别利用A、T、C、G 4 种碱基构建4条长度与原序列相等的二进制序列，用1代表与序列上相同的碱基，用0代表与序列上不同的碱基，按照A、T、C、G的顺序将4条长度与原序列相等的二进制序列连接，得到长度为原序列长度4倍的二进制序列；

②以二阶矩阵

对于上述得到的二阶矩阵，由于二阶矩阵的每一行的元素符合斐波纳契数列规则，因此，若原序列长度相等而序列不同，则该二阶矩阵一定不同。

具体地，假设序列的长度为L，每条序列由A、G、C、T这4种碱基组成，将长度为L的序列构建为二进制序列的方法为：分别按A、T、 C、G根据原序列构建4条与原序列长度相等的二进制序列，用1代表和序列上相同的碱基，用0代表和序列上不同的碱基，对于每一种碱基可以得到一条长度为L的二进制序列；然后按照A、T、C、G的顺序将4条长度为L的二进制序列连接在一起，得到1条长度为4L的二进制序列；

以序列AGTCGACAG为例，按照上述方法，其转化得到的4条二进制序列为

A：100001010

T：001000000

C：000100100

G：010010001

将其连接得到的长度为原序列长度4倍的二进制序列为 100001010001000000000100100010010001。进一步将后缀树的序列分别转化成节点数字和节点矩阵，与相应序列的长度和位置存进后缀树中，结果如图2所示。

上述方法中，如果2条序列不相同，则相应的节点数字一定不相等，基于此原则，将基于后缀树算法构建的参考基因组索引每个中间的节点序列转化成节点数字，含有“$”的序列转化成节点矩阵，实现了对参考基因组索引进行压缩的目的，显著降低了其占用内存。

以上所述的基于后缀树算法的基因组测序序列与参考基因组比对的方法中，所述将基因组测序序列与参考基因组索引进行序列比对的方法包括将节点矩阵转化为碱基序列以及将基因组测序序列与后缀树索引比对的步骤。

本发明中，节点矩阵每一行的数字符合斐波纳契数列的规则，因此，若不断将矩阵数字较大的列减去数字较小的列，重复这个步骤，最后可以得到单位矩阵

基于上述原则得到将节点矩阵转化为碱基序列的方法，该方法包括根据矩阵构建1条二进制序列，再将二进制序列转化为碱基序列的步骤。

具体地，所述二进制序列的构建包括：①首先将矩阵拆解成

②将矩阵联乘式转换为1条二进制序列：将矩阵联乘式中的矩阵

所述将二进制序列转化为碱基序列为：将新构建的二进制序列按照原序列的长度，分割成4条二进制序列，按照顺序分别对应碱基A、 T、C、G；再利用A、T、C、G 4条二进制序列构建同样长度的碱基序列，根据每个位置将对应的二进制序列的数字为1的碱基放置于对应位置，最终将节点矩阵转化成碱基序列。

具体地，以矩阵

A 10

T 00

C 00

G 01

因此，可将该节点矩阵转化为序列AG。

以上所述的方法中，将基因组测序序列与后缀树索引比对的方法包括：

①在第一层节点，查看存储节点数字的节点的原序列长度，与待比对序列长度进行比较，若待比对序列比其中一部分节点的原序列长度长，则将这部分节点原序列长度进行排序，将待比对序列首先截取最长的长度并转换成数字，与相应的节点数字进行比较，依次向下进行，直至转化的数字与某节点数字相等，将序列去掉截取部分后剩余的序列按照相同的方法与第二层节点相比；②若比较至某一层节点，待比对序列的剩余序列比该层存储节点数字的节点原序列长度都要短，则提取该层存储节点矩阵的节点，分别转换成相应碱基序列，将待比对序列与碱基序列进行比较，得到与待比对序列相同的碱基序列在节点存储的位置信息，即为待比对序列所比对到的基因组位置。

具体地，以序列“GACA”为例，将该序列往上述构建的参考基因组后缀树索引(图2)上进行比对：

①在第一层节点，序列可以拆成1个字符或7个字符，则序列 GACA只有4个字符组成，只可能拆成1个字符G，转化成数字为16.55，比对到数字16.55的节点上，剩余3个字符ACA；

②在第二个节点，有序列结束、4个字符和7个字符三种可能，则ACA可能选择4个字符和7个字符的节点，节点矩阵

另一方面，本发明还提供一种基于后缀树算法的基因组测序序列与参考基因组比对的装置，该装置包括：

获取模块，用于获取参考基因组序列；

参考基因组索引构建模块，用于根据参考基因组序列构建参考基因组索引；

序列比对模块，用于将基因组测序序列与参考基因组进行序列比对；

所述参考基因组索引构建模块包括：

后缀树构建子模块，用于构建参考基因组索引的初步后缀树；

数字转化子模块，用于将所述初步后缀树中含有分叉的节点转换成节点数字，不含有分叉的节点转换成节点矩阵，构建后续用于比对的最终后缀树。

其中，所述后缀树构建子模块具体用于：将参考基因组序列按照 100-200碱基长度、1个碱基的步长进行截断，在序列的最后插入$符号，并找到每条截断的序列的所有后缀，将各后缀按照每条后缀第一个字符在原序列中位于的坐标进行编号；将所有的后缀按照字母表的顺序进行排序；

将经过排序的每一个后缀作为节点，按照有序相邻的后缀最长公共前缀将节点合并，构成初步后缀树，每个节点记录字符串的位置。

其中，所述数字转化子模块具体用于：对于所述初步后缀树中的每条序列，分别利用A、T、C、G 4种碱基构建4条长度与原序列相等的二进制序列，用1代表与序列上相同的碱基，用0代表与序列上不同的碱基，按照A、T、C、G的顺序将4条长度与原序列相等的二进制序列连接，得到长度为原序列长度4倍的二进制序列；

以二阶矩阵

本发明提供的基于后缀树算法的基因组测序序列与参考基因组比对的装置中，所述序列比对模块包括：

碱基序列转化子模块，用于将节点矩阵转化为碱基序列；

基因组测序序列比对子模块，用于将待比对序列与后缀树索引比对。

所述碱基序列转化子模块具体用于：根据矩阵构建1条二进制序列，再将二进制序列转化为碱基序列；

所述二进制序列的构建包括：①首先将矩阵拆解成

②将矩阵联乘式转换为1条二进制序列：将矩阵联乘式中的矩阵

所述将二进制序列转化为碱基序列为：将新构建的二进制序列按照原序列的长度，分割成4条二进制序列，按照顺序分别对应碱基A、T、C、G；再利用A、T、C、G 4条二进制序列构建同样长度的碱基序列，根据每个位置将对应的二进制序列的数字为1的碱基放置于对应位置，最终将节点矩阵转化成碱基序列。

其中，所述基因组测序序列比对子模块具体用于：在第一层节点，查看存储节点数字的节点的原序列长度，与待比对序列长度进行比较，若待比对序列比其中一部分节点的原序列长度长，则将这部分节点原序列长度进行排序，将待比对序列首先截取最长的长度并转换成数字，与相应的节点数字进行比较，依次向下进行，直至转化的数字与某节点数字相等，将序列去掉截取部分后的剩余序列按照相同的方法与第二层节点相比；若比较至某一层节点，待比对序列的剩余序列比该层存储节点数字的节点原序列长度都要短，则提取该层存储节点矩阵的节点，分别转换成相应碱基序列，将待比对序列与碱基序列进行比较，得到与待比对序列相同的碱基序列在节点存储的位置信息，即为待比对序列所比对到的基因组位置。

本发明的有益效果在于：本发明提供了一种占用内存相对较小、运行速度较快的、基于后缀树算法进行序列比对的方法，有效降低了读入索引对计算机内存的要求。本发明的方法可将序列转化为对应的数字，进而将后缀树索引改为利用数字构建，通过不断将待比对序列转化成数字与索引数字进行比较，完成序列比对。

附图说明

图1为本发明实施例1中将参考基因组索引构建为后缀树的示例图。

图2为本发明实施例1中将后缀树的序列分别转化成节点数字和节点矩阵并与相应序列的长度和位置存进后缀树的示例图。

具体实施方式

以下实施例用于说明本发明，但不用来限制本发明的范围。

实施例1

本实施例提供一种基于后缀树算法的基因组测序序列与参考基因组比对的方法，利用该方法以序列“AGTCGACAG”为例构建索引后缀树，并以序列“GACA”作为待比对序列为例，将该序列往上述构建的参考基因组后缀树索引上进行比对，具体方法如下：

1、构建参考基因组索引的后缀树：

(1)在序列的最后插入“$”符号，并找出序列的所有后缀，将各后缀按照每条后缀第一个字符在原序列中位于的坐标进行编号；

1AGTCGACAG$

2GTCGACAG$

3TCGACAG$

4CGACAG$

5GACAG$

6ACAG$

7CAG$

8AG$

9G$

10$

(2)将所有的后缀按照字母表的顺序进行排序；

10$

6ACAG$

1AGTCGACAG$

8AG$

7CAG$

4CGACAG$

5GACAG$

2GTCGACAG$

9G$

3TCGACAG$

(3)构建初步后缀树；

将经过排序的每一个后缀作为节点，按照有序相邻的后缀最长公共前缀将节点合并，构成初步后缀树，每个节点记录字符串的位置(如图1所示)。

A G T C G A C A G$

1 2 3 4 5 6 7 8 9 10

2、将初步后缀树中含有分叉的节点转换成节点数字，不含有分叉的节点转换成节点矩阵，构成后续用于比对的最终后缀树：

(1)构建二进制序列

定义序列“AGTCGACAG”的长度为L，该序列由A、G、C、T 这4种碱基组成，分别按A、T、C、G根据一条序列构建4条二进制序列，用1代表和序列上相同的碱基，用0代表和序列上不同的碱基，对于每一种碱基可以得到一条长度为L的二进制序列，具体如下：

A：100001010

T：001000000

C：000100100

G：010010001

然后按照A、T、C、G的顺序将4条长度为L的二进制序列连接在一起，得到1条长度为4L的二进制序列，具体如下： 100001010001000000000100100010010001。

(2)转化成节点数字和节点矩阵

以二阶矩阵

序列“AGTCGACAG”得到的节点矩阵和节点数字如图2所示。

3、将节点矩阵转化为碱基序列

(1)根据矩阵构建1条二进制序列：

①将矩阵拆解成

②将矩阵联乘式转换为1条二进制序列：将矩阵联乘式中的矩阵

(2)将二进制序列转化为碱基序列：

将新构建的二进制序列按照原序列的长度，分割成4条二进制序列，按照顺序分别对应碱基A、T、C、G；再利用A、T、C、G 4条二进制序列构建同样长度的碱基序列，根据每个位置将对应的二进制序列的数字为1的碱基放置于对应位置，最终将节点矩阵转化成碱基序列。

按照上述方法将序列“AG”得到的节点矩阵再转化为序列，具体如下：

以矩阵

A 10

T 00

C 00

G 01

则可将节点矩阵转化为序列AG。

4、将基因组测序序列与后缀树索引比对

以序列“GACA”作为基因组测序的待比对序列，将该序列往上述构建的参考基因组后缀树索引(图2)上进行比对：

(1)在第一层节点，序列可以拆成1个字符或7个字符，则序列GACA只有4个字符组成，只可能拆成1个字符G，转化成数字为16.55，比对到数字16.55的节点上，剩余3个字符ACA；

(2)在第二个节点，有序列结束、4个字符和7个字符三种可能，则ACA可能选择4个字符和7个字符的节点，节点矩阵

由于参考基因组的序列较长，以上实施例仅以具体一条短序列作为方法操作示例，对于具体的参考基因组，可将参考基因组序列按照 150碱基长度、1个碱基的步长进行截断，在每条截断的序列的最后插入$符号，按照上述实施例的方法进行参考基因组的索引构建和基因组测序序列与参考基因组索引的比对。

虽然，上文中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。