高通量测序数据中Split read映射方法的研究

姜玥 王亚东

摘要：高通量测序技术的快速发展与广泛应用为计算机科学带来了新的挑战，read的映射问题是其中非常重要的一个部分。Split read是一类特殊的read，其出现通常是由基因组中的结构变异造成的。这类read在映射中不再保持连续序列的形式，而是包含了一定长度的空位，因此具有较高的映射难度。提出一种利用双末端测序数据的映射结果来指导split read映射的方法，这种方法可以使split read的映射难度不再与其所包含的空位数量相关，从而降低了映射过程中的搜索空间，提高映射效率。

关键词：split read； 映射； 高通量测序； 生物信息学

中图分类号：TP391 文献标识码：A文章编号：2095-2163（2013）06-0030-03

0引言

人类基因组计划的完成为人类基因组的研究提供了一套参考基因组序列，大大地简化了人类个体基因组的序列研究，因为不同人类个体基因组序列之间有着极高的相似性，现在的研究主要专注于个体基因组序列与参考基因组序列的差异，这大大地简化了研究的过程。而高通量测序技术的不断发展，则为人类基因组研究提供了有力数据支持。为了利用高通量测序数据，需要将上亿的测序短序列（read）映射到参考基因组序列上，这些read当中大部分可以以连续序列的形式被映射，但是仍有一部分read由于个体基因组序列与参考基因组序列的差异，会在映射中包含一段空位，这样的read称为split read，其映射相比于第一类read是更为困难的。Split read的映射往往可以显示个体基因组中变异区域的序列信息，对研究更快速、准确的split read映射方法有着重要的意义。

1基本概念

1.1高通量测序数据

高通量测序是一种测序DNA序列的技术。在测序过程中，将完整的样本DNA序列打碎，从中筛选出满足特定长度（通常为数百bp）的片段，然后在每个片段的一端或两端各读取一段长度为数十至数百bp的序列。这些读取出的序列长度通常远远小于被测样本DNA序列的长度，但是高通量测序技术可以同时读取大量这样的短序列，使得短序列总长度达到样本DNA长度的数倍至数十倍，从而使获得样本DNA序列成为可能。

1.2Read与split read

在高通量测序中，从打碎的DNA片段上读取出来的短序列称为read。Read是被测DNA序列的一个短片段，单个的read序列长度远远短于被测DNA序列的长度，但是通过将大量read映射到参考基因组序列的方式，就可以获得被测DNA的序列内容，如图1所示。测序时所读取的read是一段连续的序列，但是由于DNA结构变异的存在，一些read在映射结果中不再保持连续的形式，而是包含了空位，这样的read称为split read。

1.3双末端测序

在高通量测序过程中，从打碎的DNA片段的两端读取序列的方法称为双末端测序。双末端测序中获得的读取自同一片段的一对read称为一个read pair。理论上，如果被测DNA序列与参考基因组序列完全相同，read pair被映射到参考基因组之后，其中的两个read之间的距离与被测时DNA片段的长度应当是相同的。但是由于被测DNA与参考基因组序列存在差异，特别是由于结构变异的存在，read pair映射后其一对read之间的距离会与被测的DNA片段长度产生明显的差异。

2Deletion对附近read 与read pair映射所造成的影响Deletion是一种常见的结构变异形式，表现为被测DNA序列相比参考基因组序列缺失了部分序列。由于这种变异的存在，其附近的read与read pair在映射过程中会发生异常，如图2所示。从图2中可以看出，由于deletion的存在（黑色短线段），跨过deletion的read pair（左）在映射后两个read之间的距离要长于被测时两个read之间的距离，这个距离的差异恰好是deletion的长度。而跨过deletion边界的read（右）在映射时则会包含与deletion长度相同的一段空位，形成split read。

3利用read pair映射分析指导split read映射的方法目前的read映射方法出于运行效率的考虑，都会限制映射结果中所允许的空位数量与长度[1-3]。有一些利用双末端测序数据特性而特别为split read映射所设计的映射方法，利用read pair中一个映射较好的read作为基点，在临近的一段区间为另一个映射效果不好或者无法连续映射的read进行允许较多空位的映射[4]。这样的方法存在着映射效果与搜索空间相关，映射难度大，效率低等问题，如图3所示。

为了改进这些不足，本文提出一种利用deletion附近的read pair的映射结果来指导split read映射的方法。从图2中可以看出，受到deletion影响的read pair，虽然其一对read之间的映射距离发生了异常，但两个read的映射位置距离deletion的边界并不远。通过将这样存在映射异常的read pair按照映射位置与每对read之间的距离进行聚类，可以大致获得deletion边界的位置。由于split read的映射实际上只需要deletion边界处的一小段序列，而与deletion序列本身无关，因此可以每个聚类结果中的两处deletion边界位置为基点，各选择一段固定长度的序列作为参考序列进行split read映射，选择序列的长度只要确保可以包含deletion的分界点即可（图4上半部分）。通过这样的方式，split read的映射将不再与deletion本身的长度相关，因为参与split read映射的参考序列只是deletion边界处固定长度的两段序列的组合，其选取与deletion本身的长度无关。

接下来，需要将每个聚类结果附近映射效果较差或无法映射的read提取出来，这些read可能是受到了每个聚类结果所对应的deletion的影响而无法实现良好的映射，因其是候选的split read。将这些read向组合的参考序列映射需要一种序列映射算法，本文提出一种Needleman-Wunsh算法[5， 6]的变种算法来完成split read映射。变种算法同样是一种动态规划算法，其递归表达式为：

其中：

db是由两段参考基因组序列组成的横向序列，段序列的长度分别为m1和m2。qr是由read序列构成的纵向序列，长度为l。M（i，j）是当qr[i]和db[j]对齐时单元（i，j）的打分；Iqr（i，j）是qr[i]和一个空位对齐时单元（i，j）的打分；Idb（i，j）是db[j]和一个空位对齐时单元（i，j）的打分。gapopen是开始一段新空位的罚分；gapext是扩展一个空位的罚分。w（a，b）是一个打分函数，当a和b相同时打正分，反之打负分。jumpqr是matrix2中额外计算的罚分，是从matrix2中单元向matrix1中单元进行跳跃的罚分。jmax是matrix2中单元跳跃目标单元的横坐标，对于matrix2中的单元（i，j）来说，其跳跃的目标单元坐标为（i-1，jmax）。

变种算法与原算法的最大区别在于，序列比对的打分矩阵被划分为了两个部分，分别对应着deletion两个边界附近所选择出的参考序列（图4下半部分中Part 1与Part 2）。在第一部分中，全部的比对分数计算与原算法相同，在第二部分中，为每个单元计算分值时会多考虑一项，即来源于第一部分矩阵上一行中具有最高分值的单元（图4下半部分中NW-MAX单元）的打分。这个分值的计算相当于将第一部分矩阵中的部分序列比对结果与第二部分矩阵中的部分序列比对结果相连接，相连接的两个单元所在的位置就是这个映射所对应的一段连续空位的边界点。变种算法对于这种连接给出一个固定的罚分，这个罚分与两个单元的横向距离无关。在原算法中，这样的单元之间的“跳跃”是不允许的，相同的映射在原算法中需要依靠相邻单元的连续计算来完成（图4下半部分中虚线箭头所示），由于原算法中引入空位 需要罚分，因此split read的映射结果的最终分值将会受到引入的空位数量的影响，引入的空位越多，分值越低。这可能导致split read的映射结果由于引入的空位过多而导致分值过低，最终被舍弃。

4实验结果与分析

本文将所提出的算法进行程序实现，称为PRISM。通过将人类基因组中deletion注释加入到参考基因组1号染色体序列中的方式构造了一条模拟基因组序列，并使用模拟测序软件[7]对该模拟基因组序列进行模拟测序生成一套模拟数据集。在该模拟数据集上，本文将所提出的split read映射方法与一种已有的方法Pindel进行了比较。首先是运行速度上的比较，结果如表1所示。由于在取得候选split read时的标准不同，两种方法作为输入的read数量不同，但是从结果上可以看出，PRISM的输入规模略高于Pindel，而运行时间却远远短于Pindel，这证实了PRISM利用read pair分析结果来指导split read映射的方法可以大幅地提高split read映射的效率。第二项比较是split read映射效果的比较，具体结果如图5所示，可以看出PRISM在正确映射split read的能力上也要优于Pindel。

5结束语

本文提出了一种新的split read映射方法，这种方法利用split read附近的read pair映射结果分析来指导split read的映射，以达到缩小映射过程中搜索空间，提高映射效率与准确性的目的。在模拟数据实验中，通过与已有的方法进行对比，证实了本文所提出的方法在运行效率、与split read映射结果上都具有优势。

参考文献：

[1]LI H， DURBIN R. Fast and accurate short read alignment with Burrows-Wheeler transform [J]. Bioinformatics， 2009， 25（14）： 1754-1760.

[2]LANGMEAD B， SALZBERG S L. Fast gapped-read alignment with Bowtie 2 [J]. Nature methods， 2012， 9（4）： 357-359.

[3]LANGMEAD B， TRAPNELL C， POP M， et al. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome [J]. Genome biology， 2009， 10（3）： R25.

[4]YE K， SCHULZ M H， LONG Q， et al. Pindel： a pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads [J]. Bioinformatics， 2009， 25（21）： 2865-2871.

[5]DU Z H， LIN F. Improvement of the needleman-wunsch algorithm [J]. Lect Notes Artif Int， 2004， 3066：792-797.

[6]NEEDLEMAN S B， WUNSCH C D. A general method applicable to the search for similarities in the amino acid sequence of two proteins [J]. Journal of molecular biology， 1970， 48（3）： 443-453.

[7]HUANG W， LI L， MYERS J R， et al. ART： a next-generation sequencing read simulator [J]. Bioinformatics， 2012， 28（4）： 593-594.