一种基于k-mer的序列联配方法

摘要

本发明涉及一种基于k‑mer的序列联配方法。本发明通过对seq1序列和seq2序列进行k‑mer分析，获得两条序列的k‑mer集合，筛选出一致性片段。然后利用所述的一致性片段进行序列划分，进而对不同差异性片段序列进行全局联配。最后把联配的结果从5’端到3’端合并得到完整序列的联配结果。本发明利用k‑mer方法可以大大缩短序列联配时间以及联配过程中占用的计算内存。本发明建立了全新的序列联配的核心思想，并为序列联配提供了一个新的高效的技术手段。

说明书

技术领域

本发明属于生物信息领域，具体地，本发明涉及一种基于k-mer的序列比对方法。

背景技术

随着新一代测序技术454 (Roche 公司)、Solexa (Illumina 公司)和 SOLiD(ABI公司)的诞生，测序通量得到迅速提升，而测序成本急剧下降，这种突破极大地推动了基因组科学的发展。通过一代测序技术进行菌种鉴定是比传统生化鉴定更加快速、准确的鉴定方法。一代测序菌种鉴定的一般步骤就是通过检测荧光信号得到整条序列信息，然后将序列与数据库比对从而获得菌种鉴定信息。通过转录物组学和蛋白质组学等相关技术对基因表达谱、基因突变等进行匹配分析，可获得与疾病相关基因的信息。通过序列信息与基因组序列或特定基因序列(参考序列)进行联配、分析，并揭秘患病的根源。如何准确快速的从浩瀚的测序结果数据中得到基因信息的关键是序列比对。

序列比对是指通过一定算法对两条DNA或蛋白质序列进行比较，找出两者之间的最大相似性匹配。它已经成为序列比对问题和数据库搜索的基础。在现有技术中，最具有代表性的比对算法有点阵图法和动态规划算法，但是这些算法在面对大量数据时存在着处理速度慢、占用内存大等缺点。因此，对于核苷酸序列的联配方法仍需要进一步地开放和改进。

发明内容

为了克服现有技术的不足，本发明的目的是提供一种基于k-mer的序列联配方法，该方法应用于核苷酸序列联配中，能明显提高序列联配效率并减少程序运算量，快速获得两条序列的联配信息。

一种基于

步骤一，分别对输入的seq1序列和seq2序列根据第一预定长度进行

步骤二，比较步骤一所述的

步骤三，利用步骤二所述的一致性片段将序列划分成若干段差异性片段，进行差异性片段序列全局联配获得最优联配结果；如果差异性片段序列长度小于第一预定长度，则向前或向后截取一个第一预定长度的碱基并入一起联配分析；

步骤四，根据步骤三所述的最优联配结果，从5’端到3’端输出最终完整序列联配结果。

所述第一预定长度为奇数。

所述差异性片段序列全局联配包括全局比对模块和回溯模块。

所述全局比对模块实行步骤如下：

1）初始化阶段：获取待联配的subseq1序列和subseq2序列各位置上的单元信息；构建(m+1)×(n+1)的得分矩阵M，其中，m为subseq1的单元数目，n为subseq2的单元数目，subseq1序列沿顶部展开，subseq2序列沿左侧展开，得分矩阵初始化值全填充为 0；

2）计算单元得分值：用于计算得分矩阵中的单元值通过以下三个途径到达每个单元：a.来自上面的单元，代表将左侧的字符与空格比对；b.来自左侧的单元，代表将上面的字符与空格比对；c.来自左上侧的单元，代表与左侧和上面的字符比对，可能匹配也可能不匹配；即当矩阵M(i-1,j-1)、M(i,j-1)和M(i-1,j)值计算结束后，M(i,j)值才能计算；该单元的值来自于以下4个中的最大值：a. 上面单元的值-空格罚分预定值；b. 左侧单元的值-空格罚分预定值；c. 左上侧单元值+相应单元打分；d. 0；其中，所述单元值计算公式如下：

其中，

M

M

M

g表示空格罚分预定值；

E(Q

所述空格罚分预定值为-50。

所述第一预定打分值和第二预定打分值依据自定义打分矩阵。

所述单元为碱基。

所述回溯模块是根据下列步骤确定的：

(a) 确定回溯起始位置模块，所述确定回溯起始位置模块用于确定矩阵M

(b) 确定下一回溯位置模块，所述确定下一回溯位置模块用于确定基于所述回溯位置上游相邻三个位置的数值，确定下一回溯位置，其中，所述上游相邻三个位置包括行相邻位置、对角线相邻位置和列相邻位置，其中，选择数值最大的位置作为所述下一回溯位置，并且优先选择所述对角线相邻位置；如果所述最大值出现在该单元上方，则subseq1序列引入一个GAP ("-")，subseq2序列取该处碱基；如果所述最大值出现在左侧，则subseq2序列引入一个GAP ("-")，subseq1序列取该处碱基; 如果所述最大值出现在左上方，则不引入GAP，subseq1和subseq2均取相应的碱基；

(c) 重复步骤(b)，直到步骤(b)中所确定的所述下一回溯位置的行号和列号的至少之一为0；

(d) 比对结果输出模块，所述比对结果输出模块用于基于步骤(a)-(c)中所确定的回溯路线，确定所述subseq1序列与所述subseq2序列的比对结果。

所述回溯模块与所述得分矩阵单元相关联，用于基于所述得分矩阵M

其中，

M

M

M

g表示空格罚分预定值；

E(Q

所述联配结果采用三行结构输出，第一行为seq1序列，包含引入的GAP；第二行为联配结果，其中“*”表示碱基匹配，“N”表示碱基不匹配，“-”表示引入GAP；第一行为seq2序列，包含引入的GAP。

本发明的有益效果在于：

1）通过采取

2）通过全局比对，极大地满足了两条分段序列的最佳比对结果。

3）提升了序列联配的效率和准确性。

附图说明

图1为本发明实施例中的序列联配示意图。

图2 为本发明实施例中的序列分段示意图。

图3为本发明实施例中的打分矩阵和回溯模块示意图。

图4为本发明实施例中的比对结果图。

具体实施方式

为了更好的说明本发明，下面结合实施例做进一步说明，所述实施例的示例在附图中展示出。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

本发明提供的序列联配方法基于全局比对的核心思想，在此基础上进一步引入k-mer思想，依照自定义打分和特定的回溯规则，实现序列的联配。如图1所示，具体的核心步骤如下：第一步，分别对输入的seq1序列和seq2序列根据第一预定长度进行

根据本发明的实施例，上述方法进一步包括如下技术特征：

根据本发明的实施例，所述全局比对是通过如下方式进行的：获取待联配的subseq1序列和subseq2序列各位置上的单元信息；基于所述单元信息，构建(m+1)×(n+1)的得分矩阵M，其中，m为subseq1的单元数目，n为subseq2的单元数目。subseq1序列沿顶部展开，subseq2序列沿左侧展开，得分矩阵第一行和第一列数值全填充为 0；其中所述得分矩阵中的单元M

其中，M

需要说明的是，本申请所述的“比对所允许的错配数”是指在具体比对时，所允许的容错碱基数。

根据本发明的实施例，所述第一预定长度为奇数。通过模拟数据测试，发明人发现，所述第一预定长度数值并不是所有序列联配设置为统一值是最佳的，可以根据序列长度和实际情况进行选择；所述第一预定长度过长，可能会因为碱基错配导致获取不到共有k-mer，从而增加比对运算时长，比对速度较慢，也可能会影响比对结果准确性。

根据本发明的实施例，所述空格罚分预定值为-50。

根据本发明的实施例，所述第一预定打分值和第二预定打分值依据自定义 打分矩阵，如图3所示。

根据本发明的实施例，所述单元为碱基。

根据本发明的实施例，所述回溯模块（如图3所示）是根据下列步骤确定的：确定回溯起始位置模块，所述确定回溯起始位置模块用于确定矩阵M

根据本发明的实施例，所述方法可描述为：

1) 输入序列seq1和seq2进行

2) 根据

3) 根据一致性片段将序列划分成若干段差异性片段，并对分段序列进行全局比对以便获得最优联配结果。如果差异性片段序列长度小于第一预定长度，则会向前或向后截取一个

4) 根据分段最优联配结果，从5’端到3’端输出最终完整序列联配结果，如图4所示。

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