一种面向大规模基因数据的读段定位方法

摘要

本发明公开了一种面向大规模基因数据的读段定位方法；属于生物信息分析领域。该方法包括如下步骤：基因读段数据随机分割；数据的负载平衡；读段的空间索引；子读段不跨越剪切位定位；子读段跨越剪切位定位；子读段拼接；读段定位信息统计。本发明采用基于空位种子的读段定位方法，对该方法进行了改进，使之能够处理跨越剪切位读段定位的情况，设计了并行化程序并且在MapReduce框架下实现，提高了读段定位的效率。此外，本发明提出了一种负载平衡的解决方案，对个别节点执行时间较长，降低整个作业执行效率的情况具有良好的处理能力，因此具有较高的使用价值。

说明书

技术领域

本发明属于生物信息分析领域，特别涉及一种面向大规模基因数据的读段定位方法。

背景技术

基因序列匹配(读段定位)是生物信息学的一项重要研究课题，也是基因数据分析的关键 环节。基因序列需要在进行读段定位过程之后，才能展开后续一系列的生物信息分析的步骤， 如基因表达水平评估、选择性剪切事件识别、聚类等等。因此，基因读段定位已受到生物信息 领域众多研究者的广泛关注。

新一代基因测序技术的快速发展产生了海量的基因读段数据，这给传统的读段定位算法带 来了极大的挑战。庞大的基因读段数据使用传统串行算法在单核机器上运行，已经远远不能满 足研究者的需求。一般情况下定位几十GB的基因读段数据需要一周左右的时间，甚至更长， 才能获得最终的定位结果，从而导致了整个基因数据分析流程过于缓慢。目前的读段长度一般 较长，使用传统不跨越剪切位的空位种子算法可能会遗漏一些正确的定位，得到的定位结果往 往不够精确。此外，在并行化过程中总会存在某个节点的执行时间较长，使得整个作业等待个 别节点，这种情况严重影响了读段定位的效率。目前的研究热点也是针对如何解决上述问题而 展开。

典型的跨越剪切位的读段定位方法包括数据集分割、读段分割、子读段不跨越剪切位定位、 子读段跨越剪切位定位和拼接。良好的分割方法和高效的子读段跨越剪切位定位能在一定程度 上消除负载不均衡的情况以及降低算法时间复杂度，达到提高读段定位效率的效果。

发明内容

本发明为了解决基因读段定位过程中存在的效率低下的问题以及数据分割中负载不均衡 问题，提出了一种面向大规模基因数据的读段定位方法，能有效地提高基因读段定位的效率。

本发明采用的技术方案如下：

一种面向大规模基因数据的读段定位方法，包括如下步骤：

一种面向大规模基因数据的读段定位方法，其特征在于，包括如下步骤：

步骤1，基因读段数据随机分割：将给定的基因读段数据集随机分割成指定的块数；

步骤2，数据的负载平衡：利用MapReduce(映射归约)框架，从分割好的数据块中随机 抽取少量基因读段尝试执行，检测出可能执行时间较长的数据块，将其与其他数据块重新分割， 均衡负载；

步骤3，读段的空间索引：采用MapReduce框架，在Map过程中将读段分割成指定段数 的子读段，对每个子读段生成空位种子模式，将每个模式的哈希值与对应的读段信息记录在哈 希表中，构建读段的空间索引；

步骤4，子读段不跨越剪切位定位：在Map过程中利用已有的读段定位软件将子读段中能 够进行连续定位的读段首先定位，记录定位信息；

步骤5，子读段跨越剪切位定位：经过步骤4定位后，能够连续定位的子读段已经获得了 定位信息，在Map过程中根据已经定位的子读段信息，推断不能连续定位的子读段的定位信 息，实现跨越剪切位的读段定位；

步骤6，子读段拼接：经过步骤5所有读段分割成的子读段的定位信息已全部被记录，在 Map过程中将所有定位的子读段进行组装，记录所有能够连接成原始读段的定位信息；

步骤7，读段定位信息统计：在Reduce过程中将所有通过Map过程传递过来的读段定位 的信息进行汇总，统计读段定位在参考序列上的整体信息。

所述步骤2数据的负载平衡具体包括如下步骤：MapReduce框架将整个数据集默认分割为 N块，记为{S₁,S₂,...,S_N}，从每块中随机抽取一定数量的基因数据，形成相应的子块，执行读 段定位算法，检测出潜在的可能执行时间较长的分块；重新分割检测出的数据块，并且将其分 别与其他未分割的数据块合并，最后再次分割成指定的块数。

所述步骤3读段的空间索引具体包括如下步骤：在Map过程中，将对应每个数据集子块 中的读段分割为M个子读段，对每个子读段创建种子模式；每一个种子模式必然穷尽了所有 可能的含有指定错配个数碱基的方式，将每一个种子模式转换成哈希值，以键值对的形式将哈 希值和读段信息存储在哈希表里；在定位过程中就可以快速的从这张表中通过读段的哈希值获 取指定读段的信息，这张哈希表就是所有读段的空间索引。

所述步骤5子读段跨越剪切位定位具体包括如下步骤：在每一个Map中，根据读段的长 度，将子读段能够跨越的剪切位点的个数限定在指定的数量以内；根据已有的生物信息，限定 内含子最大长度，提高定位效率；根据已经定位的子读段信息，通过延展尝试匹配的方式来确 定不能够连续定位的子读段位置信息；如果原始读段的所有子读段都能够进行定位，则保留定 位信息；否则舍弃该读段。

所述步骤6子读段拼接具体包括如下步骤：首先，每个子读段都记录它能够定位到参考序 列上的所有可能位置，用一个数组存储所有位置信息；然后，从第一个子读段的第一个定位信 息，即从第一个子读段的位置数组中取出第一个位置，与第二个子读段的所有位置信息进行检 测，查看是否与其中一个位置具有顺接关系；如果有，则继续检测第二个子读段顺接的那个位 置与第三个子读段所有位置信息的关系，直至最后一个子读段；如果没有，则停止检测；接着 继续从第一个子读段的第二个位置进行检测，直至检查完第一个子读段的所有位置信息；最后， 提取所有能够完整匹配的读段定位信息。

本发明是专门针对基因读段定位而提出的方法。与现有技术相比，本发明具有以下特征：

(1)本发明针对现有基因数据量大，传统串行基因读段定位算法效率低下这一问题，提 出了并行化的读段定位算法，并且利用MapReduce框架加以实现，在不影响读段定位正确性 的条件下，提高了基因读段定位的效率；

(2)本发明对基于空位种子的不跨越剪切位的读段定位算法进行改进，使之能够处理跨 越剪切位的读段定位情况；此外，加入了生物信息：内含子的长度和剪切位点的数量，缩小了 算法的搜索空间；改进读段定位时的算法，降低了算法的时间复杂度，从而进一步提高了基因 读段定位的效率。

(3)本发明在读段定位的过程中，专门考虑了负载平衡问题，提出了一种实际的解决方 案，有效处理了个别节点执行时间长，降低整体效率的情况，因此具有较高的使用价值。

附图说明

图1为本发明整体流程图；

图2为本发明中数据的负载平衡步骤子流程图；

图3为本发明中子读段不跨越剪切位定位步骤子流程图；

图4为本发明中子读段跨越剪切位定位步骤子流程图。

具体实施方式

以下结合附图说明本发明的具体实施方式。

如图1所示，本发明公开了一种面向大规模基因数据的读段定位方法，具体步骤如下：

步骤1，基因读段数据随机分割：将给定的基因读段数据集随机分割成指定的块数；

步骤2，数据的负载平衡：采用MapReduce框架，从分割好的数据块中随机抽取少量基因 读段尝试执行，检测出可能执行时间较长的数据块，将其重新分割，分别添加到其余数据块中， 最后将其余的数据块再次进行分割，实现数据的均衡负载；

步骤3，读段的空间索引：采用MapReduce框架，在Map过程中将读段分割成指定段数 的子读段，对每个子读段生成空位种子模式，将每个模式的哈希值与对应的读段信息记录在哈 希表中，构建读段的空间索引；

步骤4，子读段不跨越剪切位定位：在Map过程中利用已有的读段定位软件将子读段中能 够进行连续定位的读段首先定位，记录定位信息；

步骤5，子读段跨越剪切位定位：在Map过程中根据已经定位的子读段信息，推断不能连 续定位的子读段的定位信息，实现跨越剪切位的读段定位；

步骤6，子读段拼接：在Map过程中将所有定位的子读段进行组装，记录所有能够连接成 原始读段的定位信息；

步骤7，读段定位信息统计：在Reduce过程中将所有通过Map过程传递过来的读段定位 的信息进行汇总，统计读段定位在参考序列上的整体信息，如定位数量、位置等；

假设读段定位过程中允许错配数量为2。

步骤1中，基因读段数据随机分割可以利用MapReduce框架实现。假设随机将读段数据 集分割为N块，则可以指定Reduce的数量为N；而在Map中可以设置一个能够产生1-N之间 数的随机产生器，每次遍历一条读段时，产生一个随机数，该随机数对应着Reduce，表示将 该条读段分配到相应的Reduce，最后只需要取出Reduce中的读段就完成数据分割了。

步骤2中，利用MapReduce框架将整个数据集默认分割为N块，记为{S₁,S₂,...,S_N}，从 每块(64MB)中随机抽取一定数量的基因数据，形成相应的子块，执行读段定位算法，检测出 潜在的可能执行时间较长的分块(假设检测出的是S₁,S₂)。

重新分割数据块，平衡负载：将S₁,S₂各自分割成N块,分别记为{S₁₁,S₁₂,...,S_1N}， {S₂₁,S₂₂,...,S_2N}。合并其余N-2块，再重新分割成N块，记为{S₃₁,S₃₂,...,S_3N}。将S₁₁,S₂₁,S₃₁合并成S′₁，按照同样的合并方法得到N个最终分块，记为{S′₁,S′₂,...,S′_N}。

其中从每块中随机抽取一定数量的样本可以通过一个MapReduce任务实现，具体步骤：

Map阶段：

(1)根据路径读取读段数据集；

(2)while(抽取的读段数量已满足要求)

①随机产生一个读段下标；

②将随机抽取到的读段送入reduce；

③删除已抽过的读段；

④抽取的读段数量加1；

(3)结束while循环。

Reduce阶段：

(1)遍历定位数量集合；

①输出参考序列名称和最终的定位数量；

(2)结束。

步骤2流程如图2所示，从各个子读段数据块中随机抽取一部分读段，做为测试数据集； 在MapReduce框架下以得到的子数据集做为输入，执行读段定位算法，以检测出可以影响整 体执行时间的个别节点；以Web浏览方式查看每一个Map节点的执行时间，根据时间的长度 推断可能负载较大的节点；将检测出的负载较大节点中的读段数据集重新分割；将重新分割的 子数据集添加到其他未分割的数据集中，达到平衡负载的目的。

步骤3中，在Map过程中，将对应每个数据集子块中的所有读段都分割为N个子读段。

针对每个子读段创建种子模式，主要思想如下：假设将子读段再分割为连续不交叉的k份， 最多允许m个错配(m<k)，对于每一个低于错配指定值的匹配，在已分的k个部分中，至少 有k-m个部分没有错配。这k-m个部分属于k个部分，因此共有组合。而可以选择的种 子模式就应该有个。每一个种子模式必然穷尽了所有可能的含有指定错配个数碱基的方 式。

将每一个种子模式转换成哈希值，以键值对的形式将哈希值和对应原始读段信息存储在哈 希表里。在定位过程中就可以快速的从这张表中通过读段的哈希值获取指定读段的信息，这张 哈希表就是所有读段的空间索引。

步骤4中，在Map过程中利用已有的读段定位软件将子读段中能够进行连续定位的读段 首先定位，记录定位信息。本发明中采用的是SeqMap软件中的读段定位算法，将读段分割为 4个子读段，错配2个，则每个子读段的种子模式就有6种。构建所有子读段的空间索引，遍 历所有参考序列，将能够连续定位的子读段定位到参考序列上。

步骤4流程如图3所示，先将子读段再次进行分割，本发明中分割为4段，创建6中空位 种子模式；然后将所有种子模式转换为哈希值，并且将哈希值与对应的读段信息存储到哈希表 中，形成所有子读段的空间索引；依次遍历参考序列中的每一个碱基；判断步骤15中形成的 哈希值在空间索引中是否存在；如果空间索引中存在该哈希值，则将对应的读段取出来与参考 序列中的碱基进行进一步匹配；如果不存在，则中断此次尝试匹配；继续遍历下个碱基。

步骤5中，经过步骤4定位后，能够连续定位的子读段已经获得了定位信息；

在每一个Map中，根据读段的长度，将子读段能够跨越的剪切位点的个数限定在指定的 数量以内(假设读段长84bp，可以限定最多允许跨越2个剪切位点)，因为跨越剪切位点数量越 多，算法就需要更多的时间去检测更多可能的定位方式；根据已有的生物信息，限定内含子最 大长度(比如已知的哺乳动物基因组中只有非常少的内含子长度大于400Kb)，这样可以减少算 法的搜索空间，提高定位效率；

根据已经定位的子读段信息，通过延展尝试匹配的方式来确定不能够连续定位的子读段位 置信息。跨越剪切位的读段定位可以分为两种情况：

第一种情况是未定位子读段的上下游子读段都已经定位，这种情况是比较容易处理的。本 发明中分别从上游读段的最后一个碱基和下游读段的第一个碱基开始匹配，分别记录下错配 0、1、2的前一个位置，在记录下的位置中查找能够连接的，且错配最少的位置，记为剪切位 点。通过这种只需一次遍历查找到剪切位点的方法，可以之间将这一过程的时间复杂度降低至 O(N)。

第二种情况是未定位子读段的上下游子段只有一个已经定位。对于这种情况，本发明中截 取未定位子读段末尾一段长度，记为h-mer；首先通过检测h-mer在上游子读段或下游子读段 中的位置，然后再利用第一种情况的处理方式进行读段定位。而在本发明中h-mer的长度取值 区别于一般算法中的取值(直接将h-mer长度固定为2)。本发明动态获取h-mer的长度，首先 通过已定位的上游或下游子读段进行延展匹配直至出现2个错配的位置，记录此时延展的长度 L，h-mer的长度就取为未定位子读段长度减去L。通过这一做法能够尽可能减少未定位子读 段尝试匹配次数，并且使得改进算法将时间复杂度O从原来的O(N³)降低至O(N²)。

进行上述跨越剪切位定位后，如果原始读段的所有子读段都能够进行定位，则保留定位信 息；否则舍弃该读段。

步骤5流程如图4所示，用已有的基于空位种子的读段定位软件将能够连续定位的子读段 定位，记录位置信息；根据已经定位的子读段信息，推断未定位的子读段；未定位子读段的可 能的第一种情况，它的上游和下游子读段均已定位成功，此时通过延展上下游子读段来确定未 定位子读段的位置信息；未定位子读段可能的第二种情况，它的上游或下游子读段之一未定位， 这种情况下可以截取未定位子读段末尾一段长度，检测它在上游子读段或下游子读段中的位 置，这样就将第二种情况转化为第一种情况了，再通过第一种情况的处理方式就可以完成定位； 拼接子读段步骤将所有子读段定位信息进行分析，拼接成完整读段定位信息。

步骤6中，经过步骤5所有读段分割成的子读段的定位信息已全部被记录；

首先，每个子读段都记录它能够定位到参考序列上的所有可能位置(可以多余1个)，可以 用一个数组存储所有位置信息；

然后，从第一个子读段的第一个定位信息，即从第一个子读段的位置数组中取出第一个位 置，与第二个子读段的所有位置信息进行检测，查看是否与其中一个位置具有顺接关系；如果 有，则继续检测第二个子读段顺接的那个位置与第三个子读段所有位置信息的关系，直至最后 一个子读段；如果没有，则停止检测；接着继续从第一个子读段的第二个位置进行检测，直至 检查完第一个子读段的所有位置信息；

最后，提取所有能够完整匹配的读段定位信息。

步骤7中，在Reduce过程中将所有通过Map过程传递过来的读段定位的信息进行汇总， 统计读段定位在参考序列上的整体信息，如定位数量、位置等。

实施例：

本实施例所使用的云计算平台共有21台服务器(1个master节点，20个slave节点)，每 个节点配置相同。节点的处理器为Intel(R)Xeon(R)CPU E5620，主频为2.40GHZ，操作系统 为64位Debian Linux操作系统。Hadoop平台版本为hadoop-0.20.2。

采用的基因数据是从拟南芥菜基因组中抽取10条基因做为参考序列，平均长度约为 12000bp。读段数据集主要有两个来源：1.用ES培养液培养拟南芥菜，重复实验3次，获得 三次情况下的拟南芥菜RNA，经测序得到3个读段样本集。2.在缺P情况下培养拟南芥菜， 同样重复实验3次，获得3个读段样本集。两种情况下获得的读段数据集共25.2GB，我们分 别称这6个读段样本集为ES1、ES2、ES3、P1、P2和P3。

本实施例包括以下部分：

1.基因读段数据随机分割：

基因读段数据随机分割可以利用MapReduce框架实现。假设随机将读段数据集分割为N 块，则可以指定Reduce的数量为N；而在Map中可以设置一个能够产生1-N之间数的随机产 生器，每次遍历一条读段时，产生一个随机数，该随机数对应着Reduce，表示将该条读段分 配到相应的Reduce，最后只需要取出Reduce中的读段就完成数据分割了。主要代码如下：

Map阶段：

(1)read_data＝get_read_dataset(read_path)；

(2)random＝Math.rand()*N+1；

(3)context.write(random,read_info)；

Reduce阶段

(1)for(count＝0；count<read_list.length；count++)

①read＝read_list.get(count)；

②context.write(task_num,read)；

(2)end.

2.数据的负载平衡：

利用MapReduce框架将整个数据集默认分割为N块，记为{S₁,S₂,...,S_N}，从每块(64MB) 中随机抽取一定数量的基因数据，形成相应的子块，执行读段定位算法，检测出潜在的可能执 行时间较长的分块。

主要代码如下：

(1)split_to_N-2(S₁,S₂,N-2)；//S₁,S₂是潜在负载较高的块；

(2)add(list_S₁,list_S₂,N-2)；//将重新分割的块添加到其余N-2块中；

(3)resplit(read_N-2_path,N)；//将N-2块重新划分为N块。

其中从每块中随机抽取一定数量的样本可以通过一个MapReduce任务实现，具体步骤：

Map阶段：

(1)read_data＝get_readdata(read_path)；

(2)while(read_num<＝fix_readNum)

①random＝Math.rand()*read_totalNum+1；

②read＝read_data.get(random)；

③context.write(fix_reduce,read)；；

④read_data.move(random)；

⑤read_num++；

(3)end while.

Reduce阶段：

(1)for(count＝0；count<read_data.length；count++)；

○1context.write(read_name,map_num)；

(2)end.

3.读段的空间索引：

本发明中读段的长度是84bp，至多允许读段跨越2个(根据生物信息学的经验在84bp长度 的读段中很少会有跨越2个剪切位点以上的读段)剪切位点。因此，在Map过程中，将对应每 个数据集子块中的所有读段都分割为4个(适中)子读段。

针对每个子读段创建种子模式，这种种子模式来自于SeqMap软件，具体可以参考文献： Jiang H,Wong W H.SeqMap:mapping massive amount of oligonucleotides to the genome[J]. Bioinformatics,2008,24(20):2395-2396.

主要代码如下：

Map阶段：

(1)list_sub_read＝split_read(read,fix_read_num)；//将读段分割为指定段数；

(2)list_seed_modes＝create_modes(list_sub_read)；//创建种子模式；

(3)for(count＝0；count<list_seed_modes.length；count++)

①hashcode＝hashcode(list_seed_modes.get(count))；

②hash_map.put(hashcode,read_info)；

//将哈希值和对应原始读段信息存储在哈希表里；

(3)end.

4.子读段不跨越剪切位定位:

在Map过程中利用已有的读段定位软件将子读段中能够进行连续定位的读段首先定位， 记录定位信息。本发明中采用的是SeqMap软件中的读段定位算法，将读段分割为4个子读段， 错配2个，则每个子读段的种子模式就有6种。构建所有子读段的空间索引，遍历所有参考序 列，将能够连续定位的子读段定位到参考序列上。

主要代码如下：

(1)for(int count＝0；count<ref_list.length；count++)//遍历参考序列数据集

①ref＝ref_list.get(count)；

②for(int loc＝0；loc<ref.length；loc++)//遍历参考序列每个碱基

a.seed_mode＝create_mode(ref.substring(loc,loc+84))；

b.if(hash_map.get(seed_mode)！＝null)match(ref,read)；

③end；

(2)end.

其中math()函数代码如下：

(1)for(int count＝0；count<read.length；count++)

①if(！read.get(count).equals(ref.get(count)))break；

else context.write(read_name,one)；

(2)end.

5.子读段跨越剪切位定位:

经过步骤4定位后，能够连续定位的子读段已经获得了定位信息；根据已经定位的子读段 信息，通过延展尝试匹配的方式来确定不能够连续定位的子读段位置信息。

主要代码如下：

(1)if(pre_subread.isFixed&&next_subread.isFixed)//第一种情况

①brige_match(read,ref)；

(2)end；

(3)if(！pre_subread.isFixed||！next_subread.isFixed))//第二种情况

①h_mer＝read.substring(read.length-L)；//动态获取h_mer

②loc＝search(h_mer,ref)；//查找h_mer位置

③brige_match(h_mer,ref)；

(4)end.

其中brige_match()函数代码如下：

(1)for(int loc＝0；lco<read.length；loc++)

①if(！read.get(loc).equals(ref.get(pre_subread.length+loc)))

a.mis_num++；

b.list_loc＝add(loc)；

②if(mis_num>2)break；

(2)end.

(3)for(int loc＝read.length；loc>＝0；loc--)

①index＝1；

②if(！read.get(loc).equals(ref.get(next_subread.loc-index)))

a.mis_num++；

b.list_loc＝add(loc)；

③if(mis_num>2)break；

(4)end,

(5)spliced_site＝find_min_mistake(list_loc).//查找错配最少的连接点

6.子读段拼接，提取所有能够完整匹配的读段定位信息。

7.统计读段定位在参考序列上的整体信息。

以上对本发明所提供的一种面向大规模基因数据的读段定位方法进行了详细介绍。值得注 意的是，具体实现该技术方案的方法和途径有很多，以上所述仅是本发明的优选实施方式，只 用于帮助理解本发明的方法及核心思想；同时，对于本领域的一般技术人员，在本发明核心思 想的基础上，做出的修改和调整都将视为本发明的保护范围。综上所述，本说明书内容不应理 解为对本发明的限制，本发明的保护范围应由所附的权利要求来限定。