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

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

Info

Publication number
CN112802553A
CN112802553A CN202011599753.0A CN202011599753A CN112802553A CN 112802553 A CN112802553 A CN 112802553A CN 202011599753 A CN202011599753 A CN 202011599753A CN 112802553 A CN112802553 A CN 112802553A
Authority
CN
China
Prior art keywords
sequence
matrix
node
length
sequences
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011599753.0A
Other languages
English (en)
Other versions
CN112802553B (zh
Inventor
张静波
姬晓勇
高司航
单光宇
伍启熹
王建伟
刘倩
唐宇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Youxun Medical Devices Co ltd
Original Assignee
Beijing Youxun Medical Devices Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Youxun Medical Devices Co ltd filed Critical Beijing Youxun Medical Devices Co ltd
Priority to CN202011599753.0A priority Critical patent/CN112802553B/zh
Publication of CN112802553A publication Critical patent/CN112802553A/zh
Application granted granted Critical
Publication of CN112802553B publication Critical patent/CN112802553B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10Sequence alignment; Homology search

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biophysics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

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

Description

一种基于后缀树算法的基因组测序序列与参考基因组比对的 方法
技术领域
本发明涉及生物信息技术领域,具体涉及一种基于后缀树算法的基因组测序序列与参考基因组比对的方法。
背景技术
序列比对是生物信息学的基本组成和重要基础。序列比对的基本思想是:基于生物学中序列决定结构、结构决定功能的普遍规律,将核酸序列和蛋白质一级结构上的序列都看成由基本字符组成的字符串,检测序列之间的相似性,发现生物序列中的功能、结构和进化的信息。
目前主要的序列比对算法分为两种:基于哈希表数据结构比对算法和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倍的二进制序列;
②以二阶矩阵
Figure BDA0002870961260000031
代表1,二阶矩阵
Figure BDA0002870961260000032
代表0,将所述长度为原序列长度4倍的二进制序列转换成矩阵联乘的运算式,并依次做矩阵的乘法,得到一个二阶矩阵,所述二阶矩阵的每一行的元素均符合斐波纳契数列规则,将所述二阶矩阵作为节点矩阵,用节点矩阵左乘权重矩阵
Figure BDA0002870961260000041
将得到的矩阵的迹作为节点数字。
对于上述得到的二阶矩阵,由于二阶矩阵的每一行的元素符合斐波纳契数列规则,因此,若原序列长度相等而序列不同,则该二阶矩阵一定不同。
具体地,假设序列的长度为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条序列不相同,则相应的节点数字一定不相等,基于此原则,将基于后缀树算法构建的参考基因组索引每个中间的节点序列转化成节点数字,含有“$”的序列转化成节点矩阵,实现了对参考基因组索引进行压缩的目的,显著降低了其占用内存。
以上所述的基于后缀树算法的基因组测序序列与参考基因组比对的方法中,所述将基因组测序序列与参考基因组索引进行序列比对的方法包括将节点矩阵转化为碱基序列以及将基因组测序序列与后缀树索引比对的步骤。
本发明中,节点矩阵每一行的数字符合斐波纳契数列的规则,因此,若不断将矩阵数字较大的列减去数字较小的列,重复这个步骤,最后可以得到单位矩阵
Figure BDA0002870961260000051
基于上述原则得到将节点矩阵转化为碱基序列的方法,该方法包括根据矩阵构建1条二进制序列,再将二进制序列转化为碱基序列的步骤。
具体地,所述二进制序列的构建包括:①首先将矩阵拆解成
Figure BDA0002870961260000052
Figure BDA0002870961260000053
联乘的形式:比较节点矩阵第一列和第二列每行对应数字的大小,若第一列对应数字比第二列对应数字大,则在矩阵联乘式开头处增加矩阵
Figure BDA0002870961260000054
并用矩阵第一列数字减去第二列数字;若第二列对应数字比第一列对应数字大,则在矩阵联乘式开头处增加矩阵
Figure BDA0002870961260000055
并用矩阵第二列数字减去第一列数字;如此反复,直至最终得到单位矩阵
Figure BDA0002870961260000056
完成将节点矩阵拆解成矩阵
Figure BDA0002870961260000057
和矩阵
Figure BDA0002870961260000058
联乘的形式;
②将矩阵联乘式转换为1条二进制序列:将矩阵联乘式中的矩阵
Figure BDA0002870961260000061
转化成数字1,矩阵
Figure BDA0002870961260000062
转化成数字0,形成二进制序列。
所述将二进制序列转化为碱基序列为:将新构建的二进制序列按照原序列的长度,分割成4条二进制序列,按照顺序分别对应碱基A、 T、C、G;再利用A、T、C、G 4条二进制序列构建同样长度的碱基序列,根据每个位置将对应的二进制序列的数字为1的碱基放置于对应位置,最终将节点矩阵转化成碱基序列。
具体地,以矩阵
Figure BDA0002870961260000063
为例,右边一列的2个数字高于左边,则用右边的数字减去左边的数字,得到矩阵
Figure BDA0002870961260000064
原矩阵与
Figure BDA0002870961260000065
右乘;第二次左边数字高于右边,用左边的数字减去右边的数字,得到矩阵
Figure BDA0002870961260000066
原矩阵与矩阵
Figure BDA0002870961260000067
右乘;以此类推,当原矩阵还原成单位矩阵
Figure BDA0002870961260000068
时,可以得到该节点矩阵是由
Figure BDA0002870961260000069
Figure BDA00028709612600000610
联乘得到,转化成二进制序列为“10000001”,由于序列由2个碱基组成,按照A、T、C、G的顺序,可以转化成:
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个字符的节点,节点矩阵
Figure BDA0002870961260000071
可转化成序列ACAG,节点矩阵
Figure BDA0002870961260000072
可转化成序列TCGACAG,则序列ACA与ACAG前三个字母相同,而与序列TCGACAG开始序列 (T)不相同,则ACA可以比对到第一个节点位置6-8上,连上前边一个碱基G,该条序列比对到了位置5-8上。
另一方面,本发明还提供一种基于后缀树算法的基因组测序序列与参考基因组比对的装置,该装置包括:
获取模块,用于获取参考基因组序列;
参考基因组索引构建模块,用于根据参考基因组序列构建参考基因组索引;
序列比对模块,用于将基因组测序序列与参考基因组进行序列比对;
所述参考基因组索引构建模块包括:
后缀树构建子模块,用于构建参考基因组索引的初步后缀树;
数字转化子模块,用于将所述初步后缀树中含有分叉的节点转换成节点数字,不含有分叉的节点转换成节点矩阵,构建后续用于比对的最终后缀树。
其中,所述后缀树构建子模块具体用于:将参考基因组序列按照 100-200碱基长度、1个碱基的步长进行截断,在序列的最后插入$符号,并找到每条截断的序列的所有后缀,将各后缀按照每条后缀第一个字符在原序列中位于的坐标进行编号;将所有的后缀按照字母表的顺序进行排序;
将经过排序的每一个后缀作为节点,按照有序相邻的后缀最长公共前缀将节点合并,构成初步后缀树,每个节点记录字符串的位置。
其中,所述数字转化子模块具体用于:对于所述初步后缀树中的每条序列,分别利用A、T、C、G 4种碱基构建4条长度与原序列相等的二进制序列,用1代表与序列上相同的碱基,用0代表与序列上不同的碱基,按照A、T、C、G的顺序将4条长度与原序列相等的二进制序列连接,得到长度为原序列长度4倍的二进制序列;
以二阶矩阵
Figure BDA0002870961260000081
代表1,二阶矩阵
Figure BDA0002870961260000082
代表0,将所述长度为原序列长度4倍的二进制序列转换成矩阵联乘的运算式,并依次做矩阵的乘法,得到一个二阶矩阵,所述二阶矩阵的每一行的元素均符合斐波纳契数列规则,将所述二阶矩阵作为节点矩阵,用节点矩阵左乘权重矩阵
Figure BDA0002870961260000091
将得到的矩阵的迹作为节点数字。
本发明提供的基于后缀树算法的基因组测序序列与参考基因组比对的装置中,所述序列比对模块包括:
碱基序列转化子模块,用于将节点矩阵转化为碱基序列;
基因组测序序列比对子模块,用于将待比对序列与后缀树索引比对。
所述碱基序列转化子模块具体用于:根据矩阵构建1条二进制序列,再将二进制序列转化为碱基序列;
所述二进制序列的构建包括:①首先将矩阵拆解成
Figure BDA0002870961260000092
Figure BDA0002870961260000093
联乘的形式:比较节点矩阵第一列和第二列每行对应数字的大小,若第一列对应数字比第二列对应数字大,则在矩阵联乘式开头处增加矩阵
Figure BDA0002870961260000094
并用矩阵第一列数字减去第二列数字;若第二列对应数字比第一列对应数字大,则在矩阵联乘式开头处增加矩阵
Figure BDA0002870961260000095
并用矩阵第二列数字减去第一列数字;如此反复,直至最终得到单位矩阵
Figure BDA0002870961260000096
完成将节点矩阵拆解成矩阵
Figure BDA0002870961260000097
和矩阵
Figure BDA0002870961260000098
联乘的形式;
②将矩阵联乘式转换为1条二进制序列:将矩阵联乘式中的矩阵
Figure BDA0002870961260000099
转化成数字1,矩阵
Figure BDA00028709612600000910
转化成数字0,形成二进制序列。
所述将二进制序列转化为碱基序列为:将新构建的二进制序列按照原序列的长度,分割成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)转化成节点数字和节点矩阵
以二阶矩阵
Figure BDA0002870961260000121
代表1,二阶矩阵
Figure BDA0002870961260000122
代表0,将步骤(1)构建的长度为原序列长度4倍的二进制序列转换成矩阵联乘的运算式,并依次做矩阵的乘法,得到一个二阶矩阵,该二阶矩阵的每一行的元素均符合斐波纳契数列规则,将该二阶矩阵作为节点矩阵;用节点矩阵左乘二阶矩阵
Figure BDA0002870961260000131
得到的矩阵的迹作为节点数字。
序列“AGTCGACAG”得到的节点矩阵和节点数字如图2所示。
3、将节点矩阵转化为碱基序列
(1)根据矩阵构建1条二进制序列:
①将矩阵拆解成
Figure BDA0002870961260000132
Figure BDA0002870961260000133
联乘的形式:比较节点矩阵第一列和第二列每行对应数字的大小,若第一列对应数字比第二列对应数字大,则在矩阵联乘式开头处增加矩阵
Figure BDA0002870961260000134
并用矩阵第一列数字减去第二列数字;若第二列对应数字比第一列对应数字大,则在矩阵联乘式开头处增加矩阵
Figure BDA0002870961260000135
并用矩阵第二列数字减去第一列数字;如此反复,直至最终得到单位矩阵
Figure BDA0002870961260000136
完成将节点矩阵拆解成矩阵
Figure BDA0002870961260000137
和矩阵
Figure BDA0002870961260000138
联乘的形式;
②将矩阵联乘式转换为1条二进制序列:将矩阵联乘式中的矩阵
Figure BDA0002870961260000139
转化成数字1,矩阵
Figure BDA00028709612600001310
转化成数字0,形成二进制序列。
(2)将二进制序列转化为碱基序列:
将新构建的二进制序列按照原序列的长度,分割成4条二进制序列,按照顺序分别对应碱基A、T、C、G;再利用A、T、C、G 4条二进制序列构建同样长度的碱基序列,根据每个位置将对应的二进制序列的数字为1的碱基放置于对应位置,最终将节点矩阵转化成碱基序列。
按照上述方法将序列“AG”得到的节点矩阵再转化为序列,具体如下:
以矩阵
Figure BDA0002870961260000141
为例。右边一列的2个数字高于左边,则用右边的数字减去左边的数字,得到矩阵
Figure BDA0002870961260000142
原矩阵与
Figure BDA0002870961260000143
右乘;第二次左边数字高于右边,用左边的数字减去右边的数字,得到矩阵
Figure BDA0002870961260000144
原矩阵与矩阵
Figure BDA0002870961260000145
右乘;以此类推,当原矩阵还原成单位矩阵
Figure BDA0002870961260000146
时,可以得到该节点矩阵是由
Figure BDA0002870961260000147
Figure BDA0002870961260000148
联乘得到,转化成二进制序列为“10000001”,由于序列由2个碱基组成,按照ATCG的顺序,可以转化成:
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个字符的节点,节点矩阵
Figure BDA0002870961260000151
可转化成序列ACAG,节点矩阵
Figure BDA0002870961260000152
可转化成序列TCGACAG,则序列ACA与ACAG前三个字母相同,而与序列TCGACAG开始序列 (T)不相同,则ACA可以比对到第一个节点位置6-8上,连上前边一个碱基G,该条序列比对到了位置5-8上。
由于参考基因组的序列较长,以上实施例仅以具体一条短序列作为方法操作示例,对于具体的参考基因组,可将参考基因组序列按照 150碱基长度、1个碱基的步长进行截断,在每条截断的序列的最后插入$符号,按照上述实施例的方法进行参考基因组的索引构建和基因组测序序列与参考基因组索引的比对。
虽然,上文中已经用一般性说明及具体实施方案对本发明作了详尽的描述,但在本发明基础上,可以对之作一些修改或改进,这对本领域技术人员而言是显而易见的。因此,在不偏离本发明精神的基础上所做的这些修改或改进,均属于本发明要求保护的范围。

Claims (10)

1.一种基于后缀树算法的基因组测序序列与参考基因组比对的方法,其特征在于,包括构建参考基因组索引以及将基因组测序序列与参考基因组索引进行序列比对的步骤,
所述构建参考基因组索引包括如下步骤:
(1)构建参考基因组索引的初步后缀树;
(2)将所述初步后缀树中含有分叉的节点转换成节点数字,不含有分叉的节点转换成节点矩阵,构建后续用于比对的最终后缀树。
2.根据权利要求1所述的方法,其特征在于,所述构建参考基因组索引的初步后缀树包括:
将参考基因组序列按照100-200碱基长度、1个碱基的步长进行截断,在每条截断的序列的最后插入$符号,并找到每条截断的序列的所有后缀,将各后缀按照每条后缀第一个字符在原序列中位于的坐标进行编号;
将所有的后缀按照字母表的顺序进行排序;
将经过排序的每一个后缀作为节点,按照有序相邻的后缀最长公共前缀将节点合并,构成初步后缀树,每个节点记录字符串的位置。
3.根据权利要求1或2所述的方法,其特征在于,所述最终后缀树的构建包括:
对于所述初步后缀树中的每条序列,分别利用A、T、C、G 4种碱基构建4条长度与原序列相等的二进制序列,用1代表与序列上相同的碱基,用0代表与序列上不同的碱基,按照A、T、C、G的顺序将4条长度与原序列相等的二进制序列连接,得到长度为原序列长度4倍的二进制序列;
以二阶矩阵
Figure FDA0002870961250000011
代表1,二阶矩阵
Figure FDA0002870961250000012
代表0,将所述长度为原序列长度4倍的二进制序列转换成矩阵联乘的运算式,并依次做矩阵的乘法,得到一个二阶矩阵,所述二阶矩阵的每一行的元素均符合斐波纳契数列规则,将所述二阶矩阵作为节点矩阵,用节点矩阵左乘权重矩阵
Figure FDA0002870961250000021
将得到的矩阵的迹作为节点数字。
4.根据权利要求1~3任一项所述的方法,其特征在于,所述将基因组测序序列与参考基因组索引进行序列比对的方法包括将节点矩阵转化为碱基序列以及将基因组测序序列与后缀树索引比对的步骤,
所述将节点矩阵转化为碱基序列的方法包括根据矩阵构建1条二进制序列,再将二进制序列转化为碱基序列的步骤;
所述二进制序列的构建包括:
将矩阵拆解成
Figure FDA0002870961250000022
Figure FDA0002870961250000023
联乘的形式:比较节点矩阵第一列和第二列每行对应数字的大小,若第一列对应数字比第二列对应数字大,则在矩阵联乘式开头处增加矩阵
Figure FDA0002870961250000024
并用矩阵第一列数字减去第二列数字;若第二列对应数字比第一列对应数字大,则在矩阵联乘式开头处增加矩阵
Figure FDA0002870961250000025
并用矩阵第二列数字减去第一列数字;如此反复,直至最终得到单位矩阵
Figure FDA0002870961250000026
完成将节点矩阵拆解成矩阵
Figure FDA0002870961250000027
和矩阵
Figure FDA0002870961250000028
联乘的形式;
将矩阵联乘式转换为1条二进制序列:将矩阵联乘式中的矩阵
Figure FDA0002870961250000029
转化成数字1,矩阵
Figure FDA00028709612500000210
转化成数字0,形成二进制序列;
所述将二进制序列转化为碱基序列为:将二进制序列按照原序列的长度,分割成4条二进制序列,按照顺序分别对应碱基A、T、C、G;再利用A、T、C、G 4条二进制序列构建同样长度的碱基序列,根据每个位置将对应的二进制序列的数字为1的碱基放置于对应位置,最终将节点矩阵转化成碱基序列。
5.根据权利要求4所述的方法,其特征在于,所述将基因组测序序列与后缀树索引比对的方法包括:
在第一层节点,查看存储节点数字的节点的原序列长度,与待比对序列长度进行比较,若待比对序列比其中一部分节点的原序列长度长,则将这部分节点原序列长度进行排序,将待比对序列首先截取最长的长度并转换成数字,与相应的节点数字进行比较,依次向下进行,直至转化的数字与某节点数字相等,将序列去掉截取部分后剩余的序列按照相同的方法与第二层节点相比;
若比较至某一层节点,待比对序列的剩余序列比该层存储节点数字的节点原序列长度都要短,则提取该层存储节点矩阵的节点,分别转换成相应碱基序列,将待比对序列与碱基序列进行比较,得到与待比对序列相同的碱基序列在节点存储的位置信息,即为待比对序列所比对到的基因组位置。
6.一种基于后缀树算法的基因组测序序列与参考基因组比对的装置,其特征在于,包括:
获取模块,用于获取参考基因组序列;
参考基因组索引构建模块,用于根据参考基因组序列构建参考基因组索引;
序列比对模块,用于将基因组测序序列与参考基因组进行序列比对;
所述参考基因组索引构建模块包括:
后缀树构建子模块,用于构建参考基因组索引的初步后缀树;
数字转化子模块,将所述初步后缀树中含有分叉的节点转换成节点数字,不含有分叉的节点转换成节点矩阵,构建后续用于比对的最终后缀树。
7.根据权利要求6所述的装置,其特征在于,所述后缀树构建子模块具体用于:将参考基因组序列按照100-200碱基长度、1个碱基的步长进行截断,在序列的最后插入$符号,并找到每条截断的序列的所有后缀,将各后缀按照每条后缀第一个字符在原序列中位于的坐标进行编号;
将所有的后缀按照字母表的顺序进行排序;
将经过排序的每一个后缀作为节点,按照有序相邻的后缀最长公共前缀将节点合并,构成初步后缀树,每个节点记录字符串的位置。
8.根据权利要求6或7所述的装置,其特征在于,所述数字转化子模块具体用于:对于所述初步后缀树中的每条序列,分别利用A、T、C、G 4种碱基构建4条长度与原序列相等的二进制序列,用1代表与序列上相同的碱基,用0代表与序列上不同的碱基,按照A、T、C、G的顺序将4条长度与原序列相等的二进制序列连接,得到长度为原序列长度4倍的二进制序列;
以二阶矩阵
Figure FDA0002870961250000041
代表1,二阶矩阵
Figure FDA0002870961250000042
代表0,将所述长度为原序列长度4倍的二进制序列转换成矩阵联乘的运算式,并依次做矩阵的乘法,得到一个二阶矩阵,所述二阶矩阵的每一行的元素均符合斐波纳契数列规则,将所述二阶矩阵作为节点矩阵,用节点矩阵左乘权重矩阵
Figure FDA0002870961250000043
将得到的矩阵的迹作为节点数字。
9.根据权利要求6~8任一项所述的装置,其特征在于,所述序列比对模块包括:
碱基序列转化子模块,用于将节点矩阵转化为碱基序列;
基因组测序序列比对子模块,用于将待比对序列与后缀树索引比对;
所述碱基序列转化子模块具体用于根据矩阵构建1条二进制序列,再将二进制序列转化为碱基序列;
所述二进制序列的构建包括:将矩阵拆解成
Figure FDA0002870961250000051
Figure FDA0002870961250000052
联乘的形式:比较节点矩阵第一列和第二列每行对应数字的大小,若第一列对应数字比第二列对应数字大,则在矩阵联乘式开头处增加矩阵
Figure FDA0002870961250000053
并用矩阵第一列数字减去第二列数字;若第二列对应数字比第一列对应数字大,则在矩阵联乘式开头处增加矩阵
Figure FDA0002870961250000054
并用矩阵第二列数字减去第一列数字;如此反复,直至最终得到单位矩阵
Figure FDA0002870961250000055
完成将节点矩阵拆解成矩阵
Figure FDA0002870961250000056
和矩阵
Figure FDA0002870961250000057
联乘的形式;
将矩阵联乘式转换为1条二进制序列:将矩阵联乘式中的矩阵
Figure FDA0002870961250000058
转化成数字1,矩阵
Figure FDA0002870961250000059
转化成数字0,形成二进制序列;
所述将二进制序列转化为碱基序列为:将二进制序列按照原序列的长度,分割成4条二进制序列,按照顺序分别对应碱基A、T、C、G;再利用A、T、C、G 4条二进制序列构建同样长度的碱基序列,根据每个位置将对应的二进制序列的数字为1的碱基放置于对应位置,最终将节点矩阵转化成碱基序列。
10.根据权利要求9所述的装置,其特征在于,所述基因组测序序列比对子模块具体用于:在第一层节点,查看存储节点数字的节点的原序列长度,与待比对序列长度进行比较,若待比对序列比其中一部分节点的原序列长度长,则将这部分节点原序列长度进行排序,将待比对序列首先截取最长的长度并转换成数字,与相应的节点数字进行比较,依次向下进行,直至转化的数字与某节点数字相等,将序列去掉截取部分后的剩余序列按照相同的方法与第二层节点相比;
若比较至某一层节点,待比对序列的剩余序列比该层存储节点数字的节点原序列长度都要短,则提取该层存储节点矩阵的节点,分别转换成相应碱基序列,将待比对序列与碱基序列进行比较,得到与待比对序列相同的碱基序列在节点存储的位置信息,即为待比对序列所比对到的基因组位置。
CN202011599753.0A 2020-12-29 2020-12-29 一种基于后缀树算法的基因组测序序列与参考基因组比对的方法 Active CN112802553B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011599753.0A CN112802553B (zh) 2020-12-29 2020-12-29 一种基于后缀树算法的基因组测序序列与参考基因组比对的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011599753.0A CN112802553B (zh) 2020-12-29 2020-12-29 一种基于后缀树算法的基因组测序序列与参考基因组比对的方法

Publications (2)

Publication Number Publication Date
CN112802553A true CN112802553A (zh) 2021-05-14
CN112802553B CN112802553B (zh) 2024-03-15

Family

ID=75805695

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011599753.0A Active CN112802553B (zh) 2020-12-29 2020-12-29 一种基于后缀树算法的基因组测序序列与参考基因组比对的方法

Country Status (1)

Country Link
CN (1) CN112802553B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114005490A (zh) * 2021-12-30 2022-02-01 北京优迅医疗器械有限公司 基于二代测序技术的循环肿瘤dna融合检测方法

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020102545A1 (en) * 1999-12-24 2002-08-01 International Business Machines Corporation, Armonk, Ny Method for changing a target array, a method for analyzing a structure, and an apparatus, a storage medium and a transmission medium therefor
CN101056993A (zh) * 2004-09-13 2007-10-17 科技研究局 用于转录作图的基因识别标签(gis)分析方法
US20130091121A1 (en) * 2011-08-09 2013-04-11 Vitaly L. GALINSKY Method for rapid assessment of similarity between sequences
EP2759952A1 (en) * 2013-01-28 2014-07-30 Hasso-Plattner-Institut für Softwaresystemtechnik GmbH Efficient genomic read alignment in an in-memory database
US20140214334A1 (en) * 2013-01-28 2014-07-31 Hasso-Plattner-Institut Fuer Softwaresystemtechnik Gmbh Efficient genomic read alignment in an in-memory database
CN104156636A (zh) * 2014-07-30 2014-11-19 中南大学 一种基于后缀数组的模糊串联重复序列识别方法
US20170109383A1 (en) * 2015-10-16 2017-04-20 Seven Bridges Genomics Inc. Biological graph or sequence serialization
CN108846016A (zh) * 2018-05-05 2018-11-20 复旦大学 一种面向中文分词的搜索算法
CN111445952A (zh) * 2020-03-25 2020-07-24 山东大学 超长基因序列的相似性快速比对方法及系统

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020102545A1 (en) * 1999-12-24 2002-08-01 International Business Machines Corporation, Armonk, Ny Method for changing a target array, a method for analyzing a structure, and an apparatus, a storage medium and a transmission medium therefor
CN101056993A (zh) * 2004-09-13 2007-10-17 科技研究局 用于转录作图的基因识别标签(gis)分析方法
US20130091121A1 (en) * 2011-08-09 2013-04-11 Vitaly L. GALINSKY Method for rapid assessment of similarity between sequences
EP2759952A1 (en) * 2013-01-28 2014-07-30 Hasso-Plattner-Institut für Softwaresystemtechnik GmbH Efficient genomic read alignment in an in-memory database
US20140214334A1 (en) * 2013-01-28 2014-07-31 Hasso-Plattner-Institut Fuer Softwaresystemtechnik Gmbh Efficient genomic read alignment in an in-memory database
CN104156636A (zh) * 2014-07-30 2014-11-19 中南大学 一种基于后缀数组的模糊串联重复序列识别方法
US20170109383A1 (en) * 2015-10-16 2017-04-20 Seven Bridges Genomics Inc. Biological graph or sequence serialization
CN108846016A (zh) * 2018-05-05 2018-11-20 复旦大学 一种面向中文分词的搜索算法
CN111445952A (zh) * 2020-03-25 2020-07-24 山东大学 超长基因序列的相似性快速比对方法及系统

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
AMOL GHOTING 等: ""Serial and parallel methods for i/o efficient suffix tree construction"", 《SIGMOD》 *
崔英博 等: ""基因组大数据变异检测算法的并行优化"", 《大数据》 *
张任文: ""生物序列索引结构的研究与实现"", 《中国优秀硕士学位论文全文数据库信息科技辑》 *
苏文鹤: ""基于后缀树的序列比对算法的设计与实现"", 《中国优秀硕士学位论文全文数据库基础科学辑》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114005490A (zh) * 2021-12-30 2022-02-01 北京优迅医疗器械有限公司 基于二代测序技术的循环肿瘤dna融合检测方法
CN114005490B (zh) * 2021-12-30 2022-04-22 北京优迅医疗器械有限公司 基于二代测序技术的循环肿瘤dna融合检测方法

Also Published As

Publication number Publication date
CN112802553B (zh) 2024-03-15

Similar Documents

Publication Publication Date Title
JP3672242B2 (ja) パターン検索方法、パターン検索装置、コンピュータプログラム及び記憶媒体
Kuruppu et al. Relative Lempel-Ziv compression of genomes for large-scale storage and retrieval
TWI636372B (zh) 用於基因定序資料的資料處理方法及系統
US11347704B2 (en) Biological graph or sequence serialization
JP5171346B2 (ja) 文字列検索システム及び方法
Haque et al. Pairwise sequence alignment algorithms: a survey
JP4912646B2 (ja) 遺伝子の転写物マッピング方法及びシステム
CN109712674B (zh) 注释数据库索引结构、快速注释遗传变异的方法及系统
JP6786144B1 (ja) Dnaに基づくデータ記憶方法、復号方法、システムと装置
US8788522B2 (en) Pair character string retrieval system
Bérard et al. Comparison of minisatellites
CN105760706A (zh) 一种二代测序数据的压缩方法
CN112802553A (zh) 一种基于后缀树算法的基因组测序序列与参考基因组比对的方法
Lin et al. Manifold de Bruijn graphs
US11594301B2 (en) DNA alignment using a hierarchical inverted index table
KR20130122816A (ko) 유전자 염기서열 압축장치 및 압축방법
CN102841988A (zh) 一种对核酸序列信息进行匹配的系统和方法
Mridula et al. Lossless segment based DNA compression
KR20190139227A (ko) K-부정합 검색을 위한 필터를 생성하는 시스템 및 방법
Dıaz-Domınguez An Index for Sequencing Reads Based on the Colored de Bruijn Graph
CN117292751A (zh) 一种基于最长路径搜索的三代序列比对方法
CN118197421A (zh) 一种针对逆补结构变异的三代序列比对方法
TW202318434A (zh) 用於處理基因定序資料的資料處理系統
Palopoli et al. Discovering frequent structured patterns from string databases: an application to biological sequences
Torkamandi Optimum Search Schemes for Approximate String Matching Using Bidirectional FM-Index

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant