CN104965999B - 一种中短基因片段测序的分析拼接方法及设备 - Google Patents

Abstract

本发明提供一种中短基因片段测序的分析拼接方法及设备，方法包括：校验read序列，去除包含错误和不可靠信息的基因序列；读取处理后的read数据，分析数据并构建k‑mer结构和四叉树结构；构建拼接存储表，记录拼接过程的进展情况和当前参与拼接的read信息；选取初始k‑mer开始拼接后，根据拼接打分公式不断选取后继k‑mer并实时更新拼接存储表结构中的信息，得到contig序列；利用read‑pair信息，基于最长公共子序列的方法合并contig序列，生成并输出super‑contig；针对以上拼接方法对性能的特殊要求，本发明设备为嵌入式手持结构；利用本发明中的方法及设备，可以快速、准确的实现中短基因片段测序的分析与拼接。

Description

一种中短基因片段测序的分析拼接方法及设备

技术领域

本发明涉及生物基因测序技术领域，特别涉及一种中短基因片段测序的分析拼接方法及设备。

背景技术

近年来，生物基因测序技术迅猛发展，虽然基因测序的精度提高，时间缩短，成本降低，但是基因测序过程中需要处理的数据量增大。因此借助计算机技术分析、处理海量数据并高效的完成基因片段拼接是基因测序的一个关键环节。

早期的Sanger第一代测序技术得到的测序片段(简称read)的主要特征是：read较长(500～1000bp)，read数量相对的较少，read之间的相互重叠关系易于发现，因此可以采用基于重叠图的贪婪算法进行拼接。对于第二代和第三代测序技术而言，read的长度只有30～50bp，read的数量更多，基于重叠图的拼接方法不再适用。目前常见的解决办法是采用基于de Bruijn图的算法，通过寻找加权图上的最优路径进行拼接。但是这种算法需要消耗大量的时间和内存来构建de Bruijn图并且容错率较低，拼接性能无法满足要求。

发明内容

为了克服上述现有技术的缺点，本发明的目的在于提供一种中短基因片段测序的分析拼接方法及设备，基于四叉树搜索方法，通过构建特殊的数据结构以及有效的拼接策略提高拼接性能，并采用最长公共子序列(LCS)方法比较两个序列的相似性，利用read-pair信息组装contig生成super-contig，本发明相应的设备可采用可视化手持嵌入式系统装置。

为了实现上述目的，本发明采用的技术方案是：

一种中短基因片段测序的分析拼接方法，包括如下步骤：

(1)、对read进行预处理操作，去掉包含错误和不可靠信息的基因序列；

(2)、读取步骤(1)中处理后的read数据，分析数据并构建k-mer结构和四叉树结构；

(3)、构建拼接存储表，记录拼接过程的进展情况和当前参与拼接的read信息；

(4)、选取初始k-mer开始拼接后，根据拼接打分公式不断选取后继k-mer并实时更新拼接存储表结构中的信息，得到contig序列；

(5)、利用read-pair信息，基于最长公共子序列的方法合并contig序列，生成并输出super-contig。

所述步骤(1)中去掉的基因序列包括：包含未能识别的碱基N的片段和A碱基含量超过90％的片段。

所述步骤(2)具体包括：

1)将read打断成长度为8的基因序列，建立k-mer数组结构；

2)初始化四叉树结构，底层采用hash表存储，通过k-mer碱基序列散列化后的值实现寻址；

3)遍历四叉树，再次读取read数据，记录每个k-mer对应read的信息。

所述步骤(2)中k-mer结构包含以下信息：

1)k-mer碱基压缩后的整数值；

2)k-mer在所有read中出现的次数；

3)k-mer在read中出现的位置以及该read的编号。

所述步骤(2)中分析数据的过程中扫描read数据两次，采用动态分配内存的方法构建k-mer结构，在第一遍扫描输入的过程中获取k-mer结构所需要存储空间的信息，在第二遍扫描输入的过程中，分配空间并且填入每个k-mer及其所在read的信息用于后续的拼接。

所述步骤(2)中构建四叉树结构时，将一个k-mer作为四叉树的父节点，四个子节点分别是后k-1个碱基补上A/C/G/T得到的k-mer。

所述步骤(3)中拼接存储表包含以下信息：

1)唯一标识一条read的字段；

2)这条read参与拼接的方向；

3)这条read参与拼接的初始位置；

4)当前k-mer出现在read中的位置；

5)拼接中断发生的位置；

6)删除标识：如果一个在决策表中的read在后继的多个碱基位上参与拼接失败，设置该标识，表示该条read参与拼接失败，当指定超时完成时，从决策表中删除设置了删除标识的read信息；

7)拼接状态标记。

所述步骤(4)中k-mer的选取遵循以下规则：

选取在所有read中出现次数最多的k-mer作为拼接的初始k-mer，当前k-mer的后k-1个碱基分别接上A/G/C/T就是四个候选后继k-mer，选取得分最高的候选后继k-mer作为当前k-mer的后继k-mer，k-mer在第i条read的得分按下式计算：

${\mathrm{score}}_i=\left\{\begin{array}{cc}k,& 1\&le;k\&le;[(L+1)/2]\\ L-k+1,& [(L+1)/2]<k<L-K+1\end{array}\right.$

其中k是k-mer在第i条read中出现的位置；L是read的长度；K是k-mer长度；[(L+1)/2]是不大于(L+1)/2的最大整数；

k-mer的总得分按下式计算：

${\mathrm{Score}}_{k-\mathrm{mer}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{score}}_i\&times;{\mathrm{num}}_i$

其中：score_i是k-mer在第i条read的得分；num_i是第i条read参与拼接的次数；n是包含该k-mer的read数量。

所述步骤(5)中利用第一代测序技术确定下来的read-pair信息获取待测片段中距离已知的两条read，当配对的两条read出现在不同的contig中时，根据read-pair信息确定contig之间的相对位置，去除重复片段合并contig，过程中采用最长公共子序列(LCS)方法比较两个序列的相似性，定义长度分别为M,N的基因序列X,Y之间的相似性函数f(X,Y)如下：

$f(X,Y)=\left\{\begin{array}{ccc}0,& \mathrm{if}& \mathrm{LLCS}(X,Y)=0\\ 1,& \mathrm{if}& \mathrm{LLCS}(X,Y)=\min (M,N)\\ \mathrm{LLCS}(X,Y)/\min (M,N),& \mathrm{other}& \end{array}\right.$

其中min(M,N)表示M,N中的较小值，LLCS(X,Y)表示X,Y间的最长公共子序列的长度，计算LLCS(X,Y)采用一种基于动态规划的递归方法，递归式如下所示：

$c[i,j]=\left\{\begin{array}{ccc}0,& \mathrm{if}& i=0\mathrm{orj}=0\\ c[i-1,j-1]+1,& \mathrm{if}& i,j>0\mathrm{and}{x}_i=y_i\\ \max \left(c\right[i-1,j],c[i,j-1\left]\right),& \mathrm{if}& i,j>0\mathrm{and}{x}_i\&NotEqual;y_i\end{array}\right.$

其中c[i,j]是序列X_i和Y_j的一个LCS的长度，X_i是基因序列X的前i个碱基组成的子序列，Y_j是基因序列Y的前j个碱基组成的子序列，通过回溯以c[i,j]为元素组成的矩阵C可以快速得到X_i和Y_j的一个LCS，矩阵C中最大的元素即最长公共子序列的长度LLCS(X,Y)。

本发明还提供了一种实现所述中短基因片段测序的分析拼接方法的设备，采用客户端-服务器硬件模型，其中本地客户端仅实现拼接过程及结果显示，数据的处理全部传输到服务器上进行，借助高性能的服务器实现基因拼接过程中数据的高速传输和海量数据的不间断处理。既可以提高拼接速度，还可以增强稳定性和可靠性，避免断电、客户端死机等异常情况对计算任务的影响。

设备硬件包括：

金属外壳1；

与金属外壳相连的液晶显示屏2；

与金属外壳相连的电源开关3、工作指示灯4和电源显示灯5；

金属外壳内部安装的实现数据快速读取的固态硬盘EMMC存储器6、双通道DDR3内存7和基于ARM Cortex-a9架构的四核处理器8；

为整个装置提供稳定电源的S5M8767电源管理单元9；

用于传输数据的USB数据模块10；

连接220V充电设备的电源输入接口11；

数据输入/输出接口12。

其中工作指示灯亮时表示设备已通电并正在工作中，电源显示灯为绿色代表电源充足，红色表示电源不足，需连接充电设备进行充电。接口11可以连接220V充电设备给本装置充电，并通过S5M8767电源管理单元9为整个装置提供稳定电源。通过输入/输出接口12传输read数据文件，经USB数据模块10传入嵌入式系统进行分析拼接处理。

与现有技术相比，利用本发明中提出的拼接方法与设备实现中短基因片段拼接测序的准确度可以达到97％，每处理1000个碱基数据仅占用1.6k内存，与现有技术相比，降低了内存消耗。

附图说明

图1是本发明中的中短基因片段测序分析拼接方法的总体流程图。

图2是本发明中的中短基因片段测序分析拼接方法的软件平台结构框图。

图3是本发明中的中短基因片段测序分析拼接装置的逻辑结构图。

具体实施方式

本发明提供了一种中短基因片段测序的分析拼接方法，其过程如图1所示。下面结合附图和实施例详细说明本发明的实施方式。

中短基因片段测序的分析拼接方法的流程图如图1所示。

1.校验read序列，将以下两种基因片段抛弃：一是包含未能识别的碱基N的片段；二是A碱基含量超过90％的片段。

2.本发明实施例中的read长度为35bp，将每条read打断成28个长度为8的k-mer(即前后两个k-mer会重叠7个碱基位的数据)，并通过构建映射关系将每个碱基压缩成一个两位的二进制数。

分析数据的过程中需要扫描read数据两次，采用动态分配内存的方法构建k-mer结构。k-mer结构需要包含的具体信息包括：k-mer碱基压缩后的整数值；k-mer在所有read中出现的次数；k-mer在read中出现的位置以及该read的编号。

在第一遍扫描输入的过程中获取k-mer结构所需要存储空间的信息，在第二遍扫描输入的过程中，分配空间并且填入每个k-mer及其所在read的信息用于后续的拼接。

构建四叉树时，将一个k-mer作为四叉树的父节点，四个子节点分别是后k-1个碱基补上A/C/G/T得到的k-mer。如:父节点是AGGTCCTA，它的四个子节点是：GGTCCTAA、GGTCCTAC，GGTCCTAG，GGTCCTAT。四叉树结构底层采用哈希表实现存储。

3.构建一种拼接存储表结构，这个结构将实时的记录拼接过程的进展情况和当前参与拼接的一些read信息，为拼接的进行和后继k-mer的选取提供参考依据。拼接存储表结构主要包括的信息有：唯一标识一条read的字段；这条read参与拼接的方向；这条read参与拼接的初始位置；当前k-mer出现在read中的位置；拼接中断发生的位置；删除标识(如果一个在拼接存储表中的read在后继的多个碱基位上参与拼接失败，设置该标识，表示该条read参与拼接失败。当指定超时完成时，从拼接存储表中删除设置了删除标识的read信息)；拼接状态标记。

4.选取在所有read中出现次数最多的k-mer作为拼接的初始k-mer，选取后继k-mer并实时更新拼接存储表结构中的信息。k-mer的选取遵循以下规则：

当前k-mer的后k-1个碱基分别接上A/G/C/T就是四个候选后继k-mer，选取得分最高的候选后继k-mer作为当前k-mer的后继k-mer。k-mer在第i条read的得分按下式计算：

${\mathrm{score}}_i=\left\{\begin{array}{cc}k,& 1\&le;k\&le;[(L+1)/2]\\ L-k+1,& [(L+1)/2]<k<L-K+1\end{array}\right.$

其中k是k-mer在第i条read中出现的位置；L是read的长度；K是k-mer长度；[(L+1)/2]是不大于[(L+1)/2]的最大整数。

k-mer的总得分按下式计算：

${\mathrm{Score}}_{k-\mathrm{mer}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{score}}_i\&times;{\mathrm{num}}_i$

其中：score_i是k-mer在第i条read的得分；num_i是第i条read参与拼接的次数；n是包含该k-mer的read数量。

按照步骤4所述方法拼接k-mer直至无法再继续拼接时，得到一条contig，重复进行步骤4所述方法可以得到多条contig。

5.利用第一代测序技术确定下来的read-pair信息获取待测片段中距离已知的两条read。当配对的两条read出现在不同的contig中时，根据read-pair信息确定contig之间的相对位置，去除重复片段合并contig，生成并输出super-contig。采用最长公共子序列(LCS)方法比较两个序列的相似性，定义长度分别为M,N的基因序列X,Y之间的相似性函数f(X,Y)如下：

$f(X,Y)=\left\{\begin{array}{ccc}0,& \mathrm{if}& \mathrm{LLCS}(X,Y)=0\\ 1,& \mathrm{if}& \mathrm{LLCS}(X,Y)=\min (M,N)\\ \mathrm{LLCS}(X,Y)/\min (M,N),& \mathrm{other}& \end{array}\right.$

其中min(M,N)表示M,N中的较小值，LLCS(X,Y)表示X,Y间的最长公共子序列的长度。计算LLCS(X,Y)采用一种基于动态规划的递归方法，递归式如下所示：

$c[i,j]=\left\{\begin{array}{ccc}0,& \mathrm{if}& i=0\mathrm{orj}=0\\ c[i-1,j-1]+1,& \mathrm{if}& i,j>0\mathrm{and}{x}_i=y_i\\ \max \left(c\right[i-1,j],c[i,j-1\left]\right),& \mathrm{if}& i,j>0\mathrm{and}{x}_i\&NotEqual;y_i\end{array}\right.$

其中c[i,j]是序列X_i和Y_j的一个LCS的长度，X_i是基因序列X的前i个碱基组成的子序列，Y_j是基因序列Y的前j个碱基组成的子序列。通过回溯以c[i,j]为元素组成的矩阵C可以快速得到X_i和Y_j的一个LCS，矩阵C中最大的元素就是最长公共子序列的长度LLCS(X,Y)。

本发明还提供了一种实现上述分析拼接方法的系统装置，包括：

在嵌入式系统中配置算法执行软件，如图2所示，软件框架包括：主窗口1、菜单区2、显示区3、参数区4、消息区5。其中菜单区包括开启拼接过程的开始6、暂停拼接过程的暂停7,、保存拼接结果的保存8和退出系统的退出9。显示区显示了当前已完成的拼接长度10和拼接完成的碱基序列11。参数区用于设置分析拼接过程中校验参数和k-mer长度。消息区5用于显示系统运行消息或提醒。

实现基因片段分析拼接方法的系统装置逻辑图如图3所示，包括：

金属外壳1、与金属外壳相连的液晶显示屏2、与金属外壳相连的电源开关3、工作指示灯4和电源显示灯5。工作指示灯亮时表示系统装置已通电并正在工作中，电源显示灯为绿色代表电源充足，红色表示电源不足，需连接充电设备进行充电。在金属外壳内部安装实现数据快速读取的固态硬盘EMMC存储器6、双通道DDR3内存7和基于ARM Cortex-a9架构的四核处理器8。接口11可以连接220V充电设备给本装置充电，并通过S5M8767电源管理单元9为整个装置提供稳定电源。通过输入/输出接口12传输read数据文件，经USB数据模块10传入嵌入式系统进行分析拼接处理。

通过本发明提出的基因片段测序分析拼接方法及装置，可以将输入长度为30～50bp的read数据集拼接至长度为10～20万bp的基因片段，并且拼接的正确率达到95％以上，从而较好的实现了中短基因片段测序的分析拼接过程。

值得注意的是，上述的具体实施方式用于解释说明本发明，仅为本发明的优选实施方案，而不是对本发明进行限制。在本发明的精神和权力要求的保护范围内，对本发明做出的任何修改、等同替换、改进等，都属于本发明的保护范围。

Claims (5)

1.一种中短基因片段测序的分析拼接方法，其特征在于，包括如下步骤：

(1)、对read进行预处理操作，去掉包含错误和不可靠信息的基因序列；

(2)、读取步骤(1)中处理后的read数据，分析数据并构建k-mer结构和四叉树结构，具体包括：

1)将read打断成长度为8的基因序列，建立k-mer数组结构；

2)初始化四叉树结构，底层采用hash表存储，通过k-mer碱基序列散列化后的值实现寻址；

3)遍历四叉树，再次读取read数据，记录每个k-mer对应read的信息；

所述k-mer结构包含以下信息：

1)k-mer碱基压缩后的整数值；

2)k-mer在所有read中出现的次数；

3)k-mer在read中出现的位置以及该read的编号；

构建四叉树结构时，将一个k-mer作为四叉树的父节点，四个子节点分别是后k-1个碱基补上A/C/G/T得到的k-mer；

(3)、构建拼接存储表，记录拼接过程的进展情况和当前参与拼接的read信息，所述拼接存储表包含以下信息：

1)唯一标识一条read的字段；

2)这条read参与拼接的方向；

3)这条read参与拼接的初始位置；

4)当前k-mer出现在read中的位置；

5)拼接中断发生的位置；

6)删除标识：如果一个在决策表中的read在后继的多个碱基位上参与拼接失败，设置该标识，表示该条read参与拼接失败，当指定超时完成时， 从决策表中删除设置了删除标识的read信息；

7)拼接状态标记；

(4)、选取初始k-mer开始拼接后，根据拼接打分公式不断选取后继k-mer并实时更新拼接存储表结构中的信息，得到contig序列；

所述k-mer的选取遵循以下规则：

选取在所有read中出现次数最多的k-mer作为拼接的初始k-mer，当前k-mer的后k-1个碱基分别接上A/G/C/T就是四个候选后继k-mer，选取得分最高的候选后继k-mer作为当前k-mer的后继k-mer，k-mer在第i条read的得分按下式计算：

其中k是k-mer在第i条read中出现的位置；L是read的长度；K是k-mer长度；[(L+1)/2]是不大于(L+1)/2的最大整数；

k-mer的总得分按下式计算：

其中：score_i是k-mer在第i条read的得分；num_i是第i条read参与拼接的次数；n是包含该k-mer的read数量；

(5)、利用read-pair信息，基于最长公共子序列的方法合并contig序列，生成并输出super-contig。

2.根据权利要求1所述中短基因片段测序的分析拼接方法，其特征在于，所述步骤(1)中去掉的基因序列包括：包含未能识别的碱基N的片段和A碱基含量超过90％的片段。

3.根据权利要求1所述中短基因片段测序的分析拼接方法，其特征在于，所述步骤(2)中分析数据的过程中扫描read数据两次，采用动态分配内存的方法构建k-mer结构，在第一遍扫描输入的过程中获取k-mer结构所需要存储空间的信息，在第二遍扫描输入的过程中，分配空间并且填入每个k-mer 及其所在read的信息用于后续的拼接。

4.一种实现权利要求1所述中短基因片段测序的分析拼接方法的设备，其特征在于：其为客户端-服务器硬件模型，其中本地客户端仅实现拼接过程及结果显示，数据的处理全部传输到服务器上进行，服务器实现基因拼接过程中数据的高速传输和海量数据的不间断处理。

5.根据权利要求4所述中短基因片段测序的分析拼接设备，其特征在于，包括：

金属外壳(1)；

与金属外壳相连的液晶显示屏(2)；

与金属外壳相连的电源开关(3)、工作指示灯(4)和电源显示灯(5)；

金属外壳内部安装的实现数据快速读取的固态硬盘EMMC存储器(6)、双通道DDR3内存(7)和基于ARM Cortex-a9架构的四核处理器(8)；

为整个装置提供稳定电源的S5M8767电源管理单元(9)；

用于传输数据的USB数据模块(10)；

连接220V充电设备的电源输入接口(11)；

数据输入/输出接口(12)。