基于MapReduce的基因组组装算法改进

何　远，张丽娜，朱兴文

（大理大学数学与计算机学院，云南大理　671003）

基于MapReduce的基因组组装算法改进

何远，张丽娜，朱兴文

（大理大学数学与计算机学院，云南大理671003）

生物信息数据的飞速增长需要新的技术引入到该学科，目前的基因组组装算法还存在着精度不高、并行化不足等缺点。对目前组装算法的分析后，提出了基于MapReduce的组装算法，通过统计去除组装过程中的错误数据，通过增加k-mer的长度消除组装过程中的重复数据，最后在MapReduce平台实现了并行组装算法，实验结果表明算法提高了组装的准确度和计算速度。

基因组组装；高通量测序；de Bruijn；MapReduce

［DOI］10. 3969 / j. issn. 2096-2266. 2016. 06. 002

生物信息学是多理论交叉的学科，它综合了计算机科学、信息技术、数学理论、统计方法等。近些年来生物数据呈指数增长，传统的分析方法已无法满足海量数据分析的要求，因此迫切需要把云计算、数据挖掘等新的计算技术应用到生物信息学领域中来〔1-2〕。基因组组装（Genome Assembly）是基因组学研究的重要内容之一，其目的是通过测序获取所研究有机体的全序列，从而实现在分子层面对有机体进行分析研究。

测序技术一般分为三代，第一代DNA测序技术〔3〕主要是指以化学降解法和双脱氧链终止法为基础的测序技术。第二代测序技术〔4〕实现高通量，大规模测序，被广泛使用。第三代测序技术〔5〕正在研究和发展的过程中，主要特点是单分子测序，最具代表性的就是Helicos公司研发的单分子测序仪和Oxford Nanopore Technologies公司正在研发的纳米单分子技术。

第三代测序技术正在研究和发展，第二代测序技术广泛应用的情况下，将云计算等新兴技术应用于生物信息学领域具有很重要的现实意义。

1　序列拼接问题研究现状

DNA测序过程大致如下〔6〕：给定一个DNA分子，将其复制n份，然后将n份DNA分子打碎，得到若干小的DNA片段。接着进行过滤，受到测序仪读取片段长度的限制，较短的片段和较长的片段被丢弃掉，这可能导致原DNA分子中有一部分区域不能被读取到，也就是说测序仪读出的序列不能完全覆盖原有的DNA分子，但只要n足够大，就可以保证原DNA分子完全被那些小片段覆盖。最后，那些留下的片段被测序仪读出，读出的DNA片段称为read。由所有read拼接成原DNA分子的过程称为基因组组装。

打碎DNA分子的过程是随机的，使得DNA片段在分子中的位置不可知，组装的过程就是重现DNA分子原有顺序的过程。组装过程中主要有序列错误、缺少足够的覆盖、重复片段等问题，重复片段是拼接过程中最棘手的问题。原因是在测序时（为保证覆盖）进行过DNA序列复制，因此无法判断它们是来自DNA序列的同一位置，还是来自不同位置的重复内容。

目前的研究主要以denove拼接为主〔7-8〕，这类算法主要包括贪心算法、Hamilton路径算法和欧拉超路算法。Hamilton路径算法没有克服重复序列的问题。欧拉超路算法由于对read进行了进一步的拆分，一定程度上克服了重复序列问题。

欧拉超路算法主要基于de Bruijn图进行拼接。但目前该算法主要是单机运算，运算速度受到限制，同时单机模式不利于完成后续的一系列信息处理。近年来云计算有了较大的发展，MapReduce技术具有计算速度更快、运算数据量更大、运算可靠性更高的特点。

通过MapReduce技术对现有组装方法改进，不仅在一定程度上克服了基因组中的重复片段造成的问题，而且利用MapReduce技术使得组装速度得到提高。

2　基因组装算法缺陷分析

序列拼接过程中的基本问题主要有：序列误差、缺少足够的覆盖、重复片段等，由测序过程可知，缺少足够的覆盖可以通过多次复制解决，序列错误往往由于实验数据错误造成，重复片段问题主要是无法区分重复片段是由于复制过程造成的，还是由于基因组自身的序列重复而造成。

考虑如下两个read〔9〕：

ATCTTATTCG

ATCTAATTCG

这两个read除了第5个位置的碱基不同外，其他部分碱基完全相同，取k-mer长度为4，其de Bruijn图如图1所示。

图1　Bubble分支结构

图1中的Bubble分支结构，主要由1位碱基错误造成，本来属于重复片段（为保证覆盖多次复制而造成），由于复制过程中有1组分子的1位碱基发生错误，使得算法误认为是新数据，从而干扰了拼接算法。

考虑下面的基因片段：

AAGACTCGACTG

取k-mer长度为4，得到的deBruijn图如图2所示。

图2　由repeat引起的分支结构

图2中由repeat引起的分支结构，其重复来自于碱基序列自身的重复，导致k-mer重复，欧拉超路算法具有一定消除重复的能力，但对repeat分支结构遍历时，将会造成拼接精度降低。图2所示的例子中，可以得到AAGACTCGACTG，或者AAGACTG，使得正确率降低。

对图2的分支结构分析后可知，其k-mer重复的原因是其长度太小（长度值为4），如果长度为5就不会出现k-mer重复。所以需要对算法进行改进。

3　MapReduce组装算法

从上述分析可知，Bubble分支结构往往由于复制过程引入单个碱基错误造成，repeat引起的分支结构往往由碱基序列自身重复引起，而碱基序列自身重复主要由所考虑的碱基数目过小造成，所以解决上述问题从以下两方面入手。

首先，对k-mer长度进行修改，其长度增加意味着重复的概率将减小，例如长度为4，则k-mer重复的概率为1/44，很容易出现重复，所以直接改为测序仪可读的最小read的碱基数，例如Solexa测序仪read长度为25～35 bp〔4〕，则直接取25为k-mer长度。则k-mer重复的概率变为1/425。

其次，由测序过程可知，不考虑碱序列自身的重复，则重复主要由复制造成，为了保证覆盖率一般复制n次，所以完成k-mer划分后，每个k-mer的数目是1-n，由概率知识可知，复制n次后，碱基片段的数目为1的概率较低，所以当k-mer数目为1时，组装过程不考虑该k-mer。但为了保证覆盖率，增加组装的准确性，如果序列出现断裂（无法完成拼接）时，使用k-mer数目为1的序列参与拼接，使其覆盖整个序列。

最后，通过MapReduce完成拼接过程的并行化，MapReduce是一种简化的分布式编程模型，核心思想是将要执行的问题分解成Map（映射）和Reduce（化简）的方式，先通过Map程序将数据切割成不相关的区块，分配（调度）给大量计算机处理，达到分布式运算的效果，再通过Reduce程序将结果汇总整理后输出。

整个拼接过程使用MapReduce改进后，具体流程如图3所示。

图3　MapReduce组装算法流程

1）将所有read进行初始化处理，read长度的最小值作为k-mer的长度，例如取25，对大于最小值的read划分为长度为最小值的k-mer，把所有k-mer平均分布到各个节点。

2）通过MapReduce并行计算，统计出k-mer的数目，统计结果大于1的数据作为下一步的输入数据，结果为1的数据保存为k-mer-num，指导后续断裂拼接。具体如图4所示。

图4　MapReduce统计k-mer原理图

3）对统计大于1的数据作为原始数据输入，对其进行编号，然后平均分派到各个节点，在Map过程中参照所有数据进行拼接，在每个节点中查找对应的后继节点。对k-mer进行de Bruijn图的遍历，算法如图5所示。头尾节点需特殊处理，图5中节点1为头节点，节点6为尾节点。

图5　MapReduce拼接流程图

4）对遍历过程中出现的断裂，通过“2）”的统计结果进行指导。如无断裂直接完成拼接。

5）最后输出对应的完整基因组序列。

4　实验及结果分析

实验平台主要由3个同配置节点利用Hadoop框架搭建的集群，每个节点拥有4 GB内存，Intel Q9400双核CPU，装载Ubuntu 12.04 LTS 64 bits系统，各节点间通过局域网互联。

4.1正确率检验BioLign〔10〕是一个针对生物信息处理的综合性软件，其中一项功能就是序列拼接，该软件的序列拼接功能采用欧拉超路算法。实验数据来自NCBI（网址是http：//www.ncbi.nlm.nih.gov）的短根瘤菌（FJ555236.1）数据，其长度为1 987 bp。

为化简实验过程，使用下载的数据对测序过程进行模拟。下载数据后，复制6次数据，使用随机函数对6组数据随机打断，丢弃小于25 bp和大于35 bp的数据，剩下的数据形成模拟基因组数据。得到模拟数据后分别使用BioLign组装方法与MapReduce算法进行组装。对组装结果和原基因序列对比，完成正确率的比较。对测序过程发生的碱基错误，通过对6组数据中的1组数据进行碱基错误模拟，然后再次组装，实验结果如表1所示。

表1　BioLign与MapReduce算法准确率比较

从结果可知，MapReduce组装方法因为去除了重复数据的影响，准确率较高。实验在理想条件下，准确率达到100%。增加k-mer的长度导致重复概率很低，以实验数据为例，其重复概率为（1987-24）/425，几乎不用考虑k-mer重复问题，当碱基数量较大时，如人类基因组达30亿bp，重复概率为（3×109-24）/425，约为1/219，从理论推导可知准确率仍达到0.999以上。

6组数据只对其中1组数据进行碱基错误模拟后能正常拼接，但拼接速度下降。因为其他5组数据的准确性不受影响，所以除了速度降低外，准确率没有下降。但BioLign组装方法因为处理Bubble分支结构和repeat引起的分支结构，其准确率有所下降，准确率主要取决于所采用的策略。

4.2拼接速度检验对单机与MapReduce组装方法的运行速度进行实验，实验对数据初始化、统计k-mer部分的操作时间等统一没有计算。实验仍以模拟基因组数据进行；对测序过程发生的碱基错误，仍通过对6组数据中的1组数据进行碱基错误模拟，获得数据后分别进行组装，实验结果如表2所示。

表2　单机与MapReduce组装方法速度比较结果

从结果可看出MapReduce组装方法的速度接近单机的3倍，对错误数据处理，因为需要满足覆盖率的要求，对断裂部分特殊处理，增加了运行时间，但速度仍接近3倍。当数据错误率增加时，MapReduce组装方法的耗时会增加。同时，在拼接过程中，可能会因为节点数据的差异造成负载不均衡问题而导致拼接速度的下降。

5　结束语

MapReduce算法具有更高的准确性，更快的拼接速度。但对此方法仍有一定不足，如实验过程使用了理想化数据。组装过程需要对数据进行预处理，中间需要通过统计结果指导断裂处的拼接，如果断裂较多会造成多次拼接操作，会严重影响拼接速度。此外需要对头尾节点进行特殊处理等工作消耗的时间也不少，如何更高效的提高精度和速度需要进一步研究。

〔1〕杨帅，胡宗倩，伯晓晨，等.云计算在生物医学中的应用〔J〕.中国科学，2013，43（7）：569-578.

〔2〕SCHADT E E，LINDERMAN M D，SORENSON J，et al. Computational solutions to large-scale data management and analysis〔J〕.Nat Rev Genet，2010，11（11）：647-657.

〔3〕MAXAM A M，GILBERT W.A new method for sequencing DNA〔J〕.Proc Natl Acad Sci USA，1977，74（3）：560-564.

〔4〕MILNE I，BAYER M，CARDIE L，et al.Tablet-next generation sequence assembly visualization〔J〕.Bioinformations，2010，26（3）：401-402.

〔5〕CLARKE J，WU H C，JAYASINGHE L，et al.Continuous base identification for single-molecule nanopore DNA sequencing〔J〕.Nat Nanotechnol，2009，4（4）：265-270.

〔6〕王旭.基于de Brujin图的DNA Contig生成算法〔D〕.哈尔滨：哈尔滨工业大学，2011：6.

〔7〕孙晓斐.基因组序列 de novo拼接系统的设计与实现〔D〕.哈尔滨：哈尔滨工业大学，2014：5.

〔8〕王小艳.基于De Bruijn图的基因拼接算法研究〔D〕.武汉：武汉理工大学，2014：5.

〔9〕王东阳，任世军，王亚东.DNA序列拼接中de Bruijn图结构的研究〔J〕.智能计算机与应用，2011，8（4）：20-27.

〔10〕FERRARINI M，MORETTO M，WARD J A，et al.An evaluation of the PacBio RS platform for sequencing and de novo assembly of a chloroplast genome〔J〕.BMC Genom，2013，14（1）：670-671.

〔Abstract〕The rapid increase of biological information data requires the import of new technology.At present，the genome assembly algorithm is neither precised nor parallelize.A new algorithm based on MapReduce is proposed after analysis of the current assembly algorithm.The error data is removed through statistics way，and the duplicate data is eliminated by increasing the length of the k-mer in the process of assembly.Finally，the parallel assembly algorithm is realized in MapReduce platform.The experimental results show that the accuracy and speed of this algorithm are improved.

〔Key words〕genome assembly；high throughput sequencing；de Bruijn；MapReduce

（责任编辑袁霞）

Improvement of Genome Assembly Algorithm Based on MapReduce

He Yuan，Zhang Lina，Zhu Xingwen
（College of Mathematics and Computer，Dali University，Dali，Yunnan 671003，China）

TP309.2：TP399

A

2096-2266（2016）06-0004-04

大理大学青年教师科研基金资助项目（KYQN201218）

2016-01-25

2016-03-02

何远，实验师，主要从事计算机科学研究.