英语原文共 7 页，剩余内容已隐藏，支付完成后下载完整资料
用Burrows–Wheeler变换进行快速准确的短序列比对
Heng Li and Richard Durbin
Wellcome Trust Sanger Institute, Wellcome Trust Genome Campus, Cambridge, CB10 1SA, UK
Received on February 20, 2009; revised on May 6, 2009; accepted on May 12, 2009
Advance Access publication May 18, 2009
Associate Editor: John Quackenbush
摘要
动机：新DNA测序技术产生了大量短序列，需要开发快速准确的序列比对程序来处理这些数据。第一代基于哈希表格法的程序已开发完成，它包括MAQ，MAQ准确、功能丰富并且足够快来比对单个体的短序列。然而，MAQ不支持单端序列的缺口比对，这使得它不适合较长序列的比对，此处经常发生插入缺失当比对数百个个体的基因测序时，MAQ的速度就是个问题了。
结果：我们实现了一个新的序列比对安装包——Burrows-Wheeler比对工具（BWA），它基于逆向搜索通过Burrows–Wheeler变换（BWT）来有效的比对一组大的参考序列（如人类基因组）中的短序列（包含错配和缺失）BWA支持基础空间序列（如通过Illumina测序机）和色彩空间序列（如通过AB固体机）仿真和实际数据的评估显示BWA比MAQ快sim;10–20times;且准确性差不多。此外，BWA输出对比为新标准SAM（序列对比/图）格式。比对后的变异调用和其他逆向分析可被实现通过开源SAM工具软件包。
网站: http://maq.sourceforge.net
联系我们: rd@sanger.ac.uk
1 介绍
Illumina/Solexa测序技术在一次机器运行中通常产生50–200百万 32–100 bp的序列。比对大量短序列到一个和人类基因组等同大小的基因组，对现有序列比对程序构成了一个巨大挑战。为了满足高效、准确的短序列映射，开发出了许多新的比对程序。其中有Eland, RMAP, MAQ, ZOOM, SeqMap, CloudBurst和SHRiMP，它们通过散列短序列并对参考序列进行扫描，这类中的程序通常具有灵活的内存占用，但在比对很少的reads时可能需要从头开始扫描整个基因组。第二类软件包括SOAPv1 (Li et al. 2008b), PASS (Campagna et al., 2009), MOM (Eaves and Gao, 2009), ProbeMatch (Jung Kim et al., 2009), NovoAlign (http://www.novocraft.com), ReSEQ (http://code.google.com/p/ re-seq), Mosaik (http://bioinformatics.bc.edu/marthlab/Mosaik) 和 BFAST (http://genome.ucla.edu/bfast), 散列基因组，这类程序可以使用多线程进行并行处理，但它们通常需要大量的内存来构建人类基因组的索引。此外，这些软件通常采用迭代法，这使它们的速度受限于测序错误率。第三类包括slider (Malhis et al., 2009)，它进行对比通过归并排序参考子序列和短序列。
近日，采用Burrows–Wheeler变换（BWT）的基于字符串匹配的理论引起了一些机构的兴趣，这些机构是SOAPv2，Bowtie 和 BWA（我们这篇文章中介绍的新对比方式）的开发者，本质上，用BWT使用逆向搜索，我们可以有效模仿基因组自上而下全部的前缀特里结构且使用相对较小的内存占用，不管基因组大小如何，可在O(m)时间内测量精确匹配长度为m的基因串，为了模糊搜索， BWA从隐式前缀特里结构中抽取出与被查询的短序列编辑距离小于k的明确子串，由于前缀特里结构一条路径确切的重复会被压缩，我们不需要对比短序列和每个重复副本，这是基于BWT的算法这么有效的主要原因。
在本文中，我们将对BWT的背景和精确匹配时的逆向搜索给予充分的介绍，并提出BWA中使用的模糊匹配算法。我们通过对比模拟数据和真实数据、检查成对一致的序列映射成分以及通过计算一个混合基因组中的未比对序列映射得到真实末端匹配数据来评估BWA的表现
2 方法
2.1 前缀特里结构和字符串匹配
X字符串的前缀特里结构是树状的，它的每个边都标有一个符号，并拼接在从叶到根的路径的边上的符号的并置字符串给出X的唯一前缀。在前缀特里结构中，从节点到根的边缘符号的并置字符串给出了唯一X子串，它由节点表示为字符串。需要注意的是X的前缀特里结构同反向X的后缀特里结构是相同的，因此后缀特里结构理论同样适用于前缀特里结构。通过前缀特里结构，测试一个查询W是否是X的确切子串相当于寻找代表W的节点，它能在O(|W|)时间内通过比对各个W中从根开始到边缘的符号完成。为了允许不匹配情况，我们可以彻底遍历特里结构并比对每个可能路径的W。稍后我们将介绍如何使用W的前缀信息来加速搜索。图一给出了一个lsquo;GOOGOLrsquo;前缀特里结构的例子。每个节点中的后缀列阵（SA）间隔在2.3节进行说明。
图 1. 在一个节点中的两个数字表示由节点所表示的字符串的SA的间隔（参见第2.3节）虚线显示了查询字符串lsquo;LOLrsquo;强力搜索的路径，至多允许一个错配。方形标注的边上的标签指搜索中不匹配的查询。唯一命中的是加粗节点1.1，它代表字符串lsquo;GOLrsquo;
2.2 Burrows–Wheeler 变换
假设！是一个字母表，符号$不在字母表！中，而且按照字典序排列，$比字母表！中的任何字母都要小。任何字符串X= a0a1 ...anminus;1都以符号$结尾（即anminus;1=$）且$仅在字符串结尾出现。X[i]=ai , i=0,1,...,nminus;1 表示字符串X的第i个字符，X[i,j]=ai ...aj 表示字符串X的一个子串，Xi=X[i,nminus;1] 表示字符串X的一个后缀。X的后缀数组S为整数0...nminus;1的排列，比如S(i)就代表第i个最小后缀的初始位置。字符串X的BW变换定义为：当S(i)=0时，B[i]=$，S(i)0时，B[i]=X[S(i)minus;1]。|X|表示字符串X的长度，所以|X|= |B|=n。构建BWT和后缀列阵的方法如图二所示。
此算法所需的时间和空间都与问题规模成平方关系。不过这不是必须的，实际上我们通常先构造一个后缀数组再进行BW变换。大部分算法构造后缀数组至少需要nlceil;log2nrceil;位的工作空间，对于人类基因组来说大概需要12GB。Hon et al. (2007) 提出了一个新的算法，只用了n位的工作空间而且对于人类基因组的BW变换峰值内存消耗lt;1GB。BWT-SW (Lam et al., 2008) 实现了这个算法，我们修改了它的源代码使其适用于BWA。
图 2. 构造字后缀数组以及对字符串X=googol$进行BW变换。字符串X循环移位得到7个按字典序排序的字符串，排序后得到后缀数组的第一个字符的位置(6,3,0,5,2,4,1)和循环得到的字符串的最后一位的串联得到BWT字符串lo$oogg.
2.3 后缀数组区间和序列比对
如果字符串W是X的子串，那么W在X中的每一个出现的位置都会在X排序后的后缀数组的一个区间中，这是因为所有以W为前缀的X的后缀都被排序到了一起，更具这个现象，我们进行如下定义：
在特殊的情况中，W是一个空的字符串，此时。区间称为W的SA区间，W在X上每一个出现的位置的集合是。在图2的例子中，字符串lsquo;gorsquo; 的SA区间是[1,2] 后缀数组在这个区间的值是3和0表示的是字符串lsquo;gorsquo;所有的出现的位置，当知道了后缀数组的区间的时候，我们就能得到位置。因此，序列比对等价于查找X的子串的SA区间。对于精确匹配，我们能找出唯一的区间，而模糊匹配会得到很多个区间。
2.4 精确匹配：反向搜索
aisin;!，C(a)表示X[0,nminus;2]中字典序小于a的字符个数。O(a, i)表示字符a在B[0,i]中出现的次数。Ferragina 和 Manzini (2000)证明了如果W是X的子串，那么：
当且仅当aW是X的子串时，。这个结论使得可以在O(|W|)时间内迭代计算W以及W的后缀的R来得到W在X上出现的次数，这个过程称为反向搜索。值得注意的是等式(3)(4)实际上实现了对于给定的前缀树X，如果我们知道父节点的SA区间，就可以通过自上而下遍历前缀树在固定的时间内计算出子节点的SA区间。从这个层面上来说，反向搜索等价于前缀树的精确匹配，但是没有显式将前缀树保存在内存中。
图3. 模糊搜索子串W的SA区间的算法。参考字符串X以$结尾，W以A/C/G/T结尾。函数InexactSearch(W,z)返回差异（错配或缺失）个数不超过z的子串W的SA区间。InexRecur(W,i,z,k,l)在后缀Wi 1在区间[k, l]上时，递归的计算差异个数不超过z的子串W[0,i]的SA区间。以星号开头的两行分别是从X插入和删除。D(i)是字符串W[0,i]的差异个数的下界。
2.5 模糊匹配：有界遍历/回溯
图3 给出了一个递归算法来搜索差异（错配或缺失）个数不超过z的X的子串W。本质上，这个算法使用反向搜索从基因组中取出不同的子串。这个过程是由数组D(·)界定的，其中D(i)是字符串W[0,i]的差异个数的下界。对D估计的越好，算法搜索的空间就越小，效率也越高。对于所有的i，令D(i)=0来可以到一个自然的边界，但是这个算法效率很低而且结果的差异个数显然是指数增长。
图4. 计算D(i)的等价算法
图3所示的函数CalculateD提供了一个虽然不是最好的，但是也比较好的边界。它在概念上等同于与更容易理解的图4所示的函数CalculateD。我们使用未补全的反向参考序列的BW变换来检测W的子串是否也是X的一个子串。需要注意的是单独测试BW变换字符串B的CalculateD的复杂度是O(|W|2)而不是图三所示的O(|W|)。
为了说明D的作用，我们再回过头来看在图一所示的X=GOOGOL$例子中搜索字符串W =LOL，如果我们使对于所有的i，D(i)=0，且禁止跳空匹配（删除星号开头的两行），InexRecur的调用图是一个树状的，图1中的虚线指示了搜索的路线，然而对于CalculateD，我们知道当D(0)=0且D(1)=D(2)=1时，我们可以避免陷入G和O的子前缀树从而得到小得多的搜索空间。
图3所示的算法在理论上可以保证找差异个数不超过z的所有区间，在实践中我们还进行了一些修改：首先我们制造了一些gap opens以及gap extensions，使得更加接近于真实的生物数据，其次我们使用了堆状的数据结构，而不是递归。堆状结构被优先在局部命中
剩余内容已隐藏，支付完成后下载完整资料
资料编号：[148389]，资料为PDF文档或Word文档，PDF文档可免费转换为Word
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
