快速检测复杂变体和短阅读中剪接的SNP耐受性

摘要

动机：下一代测序捕捉了相对于参考基因组或转录组的读数的序列差异，包括剪接事件和涉及多个错配和长INDEL的复杂变体。基于从基因组索引中合并和筛选位置列表的连续约束搜索过程，我们提出了快速检测复杂变体和短阅读中剪接的计算方法。我们的方法是在GSNAP（Genomic Short-read Nucleotide Alignment Program）中实现的，该程序可以对齐短至14nt、任意长度的单端和双端读取。它可以使用概率模型或已知剪接位点的数据库检测单个读取中的短剪接和长剪接，包括染色体间剪接。我们的程序还允许对所有可能的主要和次要等位基因组合的参考空间进行SNP耐受性比对，并可以对亚硫酸氢盐处理DNA的读取进行比对，以研究甲基化状态。

结果：在比较测试中，GSNAP的速度与现有程序相当，尤其是在读取≥70 nt时，并且在检测具有四个或更多不匹配或插入1-9 nt和删除1-30 nt的复杂变体时最快。尽管SNP容差并没有显著增加对齐率，它影响7-8%的转录读取的比对结果，通常是通过显示读取的备用基因组映射。亚硫酸氢盐转化DNA的模拟显示，在36 nt读取中，6%的读取和70 nt读取中3%的读取中，唯一识别基因组位置的能力下降。

可利用性：C语言中的源代码和Perl语言中的实用程序可以作为GMAP包的一部分免费下载，网址为http://share.gene.com/gmap.

联系人： twu@gene.com

1简介

迄今为止，已经开发了许多程序，用于将下一代测序技术的短读数与参考基因组或转录组进行比对，如Illumina/Selexa（美国加利福尼亚州海沃德）和ABI/SOLiD（美国加州福斯特城）。由于给定的样本可以产生大量的短读，因此重点放在了速度上。因此，最近的项目，如Bowtie（Langmead等。,2009)、BWA（Li和Durbin，2009)和SOAP2（Li等。,2009)，展示了后缀数组（Manber和Myers，1993)，使用Burrows–Wheeler变换（BWT；BurrowsandWheeler，1994)，可以快速映射精确匹配的读取，或者相对于引用有一些不匹配或短插入或删除（索引）。

除了速度之外，扩大在阅读中可以检测到的可能变体的范围也很重要，因为有趣的生物学可能不仅表现为单核苷酸多态性（SNP）或参考突变，还表现为更复杂的现象，例如多重错配，长indels及其组合。这种复杂的变异是遗传多样性的重要来源。例如，indels代表7-8%的人类多态性，25%的编码indels长于3nt（Bhangale等。,2005; 韦伯等。,2002). 影响多种氨基酸的长indels可能会产生重大的生物后果。此外，随着阅读继续延长，从最初的～30 nt到目前的75–100 nt，阅读更有可能与参考文献存在多重或复杂差异，这使得复杂变体的检测变得更加关键。

剪接事件还产生了其他重要的生物学现象，这些现象提供了对基因结构、选择性剪接、基因融合和染色体重排的见解。虽然剪接可以使用通用的基因组定位和比对程序（如BLAT（Kent，2002)或GMAP（Wu和Watanabe，2005)短阅读带来了挑战，因为它们经常与基因组中的许多位置对齐，并且因为它们通常在外显子-外显子连接的一端或两端缺乏足够的序列特异性，无法准确定义连接。

检测短阅读中剪接的一种解决方案是将它们与参考转录组对齐，可能会用人工构建的外显子-外显子片段进行扩增（Wang等。,2008). 然而，这种方法将只识别已知或预测的外显子组合，而不是通过外显子跳过、隐式剪接或基因融合发生的意外外显子对。TopHat项目采用的另一种方法（Trapnell等。,2009)，分析映射读取的整个数据集，以识别给定邻域中外显子之间的剪接位点连接。然而，这种方法要求外显子具有足够高的表达，并将错过由低水平的个别阅读跨越的剪接事件。QPALMA项目（博纳）提供的第三种方法等。,2008)，可以使用Smith–Waterman型比对和经过专门训练的剪接位点模型，对齐外显子-外显子连接处的单个读取。所有这些方法都局限于仅识别局部外显子-外显子连接，而不是意外的远距离或染色体间基因融合事件。

为了扩大可以从短阅读中推断出的生物现象的范围，我们开发了快速且记忆有效的方法来检测单个阅读中的复杂变体和剪接。我们的方法在GSNAP公司（基因组短阅读核苷酸校准计划），它可以校准单个和成对阅读，短至14nt，长至所需长度。我们的程序可以检测到拼接、多个不匹配、长索引及其组合，最多可检测到用户指定的总点数，每次读取仅限于一个拼接或索引，前提是读取（或索引或拼接每一端的部分读取）具有与参考序列匹配的连续14 nt延伸(图1A） ●●●●。（该程序的未来版本可能允许每次读取多个索引或可能多个拼接。）

我们的程序可以使用两种证据识别短读中的拼接。首先，它可以使用供体和受体剪接位点的概率模型来评估周围的基因组序列(图1B） ●●●●。其次，它可以利用用户提供的已知外显子-内含子边界的数据库，避免了概率模型的假阳性和假阴性结果(图1C） ●●●●。已知的剪接位点将揭示大多数选择性剪接和基因融合事件，这些事件通常通过外显子跳跃和内含子之间的交叉发生。GSNAP可以依靠这两种证据中的一种或两种来识别剪接事件，包括那些不匹配的剪接事件以及部分剪接或“半内含子”事件，其中一端有足够的序列与一个外显子对齐，而另一端缺乏足够的序列来识别另一个外隐子(图1C） ●●●●。剪接事件不仅可以跨越正常剪接或选择性剪接中的短距离，还可以跨越长距离染色体内缺失或反转以及染色体间易位(图1D） ●●●●。

我们的程序还实现了将读取内容与单个参考序列对齐的功能，而且还实现了与数据库（如dbSNP）中所有可能的主要和次要等位基因组合的参考“空间”对齐的功能（Sherry等。,2001). 通过对齐参考空间而不是单个参考序列，我们的程序避免了将次要等位基因视为不匹配，从而在对齐过程中惩罚这些基因型。一个例子说明了耐SNP比对的实用性，该例子最初建议一个样本特有的剪接连接，但实际上是由于附近的两个次要等位基因导致读取无法对齐(图1E） ●●●●。在耐SNP比对中，次要等位基因被视为与参考空间的匹配，而不是与参考序列的不匹配。这个想法已经在其他地方得到了实施（曼斯克和奎亚特科夫斯基，2009)，但在某种程度上需要29GB的内存才能与人类基因组进行SNP耐受性比对。在本文中，我们展示了如何使用更小的内存需求来执行此功能。

我们的方法可以推广到其他任务中，例如从用亚硫酸氢钠（BS）处理的DNA中绘制读取数据，以研究甲基化状态（Lister和Ecker，2009). BS将基因组DNA中的每一个非甲基化胞嘧啶转化为尿嘧啶，尿嘧啶随后以胸腺嘧啶的形式出现。从BS-converted DNA读取的错误率很高，这为有效检测多个不匹配提供了额外的动机。尽管现有的比对程序可以用来处理这些数据，但通过将参考序列和读数中的所有胞嘧啶转换为胸腺嘧啶（Deng等。,2009)一个微妙的问题是，这种方法掩盖了参考胸腺嘧啶和阅读胞嘧啶之间的不匹配。GSNAP中的数据结构使其能够对照参考序列或SNP耐受参考空间，将BS-seq读数与基因组T至读数-C错配的显式检测进行比对。

2方法

2.1概述

我们将比对视为参考序列中基因组区域空间上的搜索问题，或者如果允许空白，则将其视为区域组合。（虽然参考序列可能由染色体、连接体、转录本或人工片段组成，但我们通过将其称为“基因组”来简化我们的论述。）搜索包括生成、筛选和验证候选基因组区域的步骤，其效率取决于设计生成和过滤步骤，以产生尽可能少的候选。包括MAQ（Li等。,2008年a)、RMAP（史密斯等。,2008)，SeqMap（Jiang和Wong，2008)和RazerS（Weese等。,2009)，对读取进行预处理，然后根据基因组扫描读取索引，生成并筛选候选基因组区域。

对于大型基因组，预处理基因组比读取创建基因组索引文件更有效，后者为给定的寡聚体提供基因组位置。基因组索引还允许部分读取与任意基因组区域对齐，这是远距离剪接检测所需的。每个引用序列只需进行一次索引，生成的索引文件可供每个新数据集使用。所有长度的低聚物都可以使用后缀数组或其压缩的BWT等价物进行索引，如Bowtie、BWA和SOAP2中使用的那样，它们可以紧凑地表示一个参考序列，对于人类大小的30亿nt基因组，以2GB表示。

然而，当只有单一低聚物长度q个算法需要一个简单的哈希表（Ning等。,2001)或q个-克指数（拉斯穆森等。,2006)应用于基因组就足够了(图2A） ●●●●。此数据结构由所有可能的偏移文件（或查找表）组成q个-mers，带有指向位置文件（或出现表）的指针，其中包含每个位置的基因组位置列表q个-梅尔。为了使我们的搜索算法最有效地工作，重要的是位置文件中的每个位置列表都要预先分类，这样可以快速计算多个位置之间的交点q个-mer查找。交集过程还要求每个位置列表中的位置在运行时根据其在给定读取中的位置进行调整，因此它们对应于基因组和读取之间对齐矩阵中的对角线。虽然我们的比对算法可能与另一个数据结构一起工作，该数据结构为给定的q个-mer，后缀数组需要在运行时对每个位置列表进行排序。

一套n个排序后的列表可以及时合并O（运行）(我日志n个)，其中我是列表长度的总和，使用基于堆的多路合并程序（Knuth，1973). 合并过程不仅可以产生候选基因组区域的列表，还可以产生关于支持每个区域的位置列表的计数和读取位置的信息。此支持信息可以提供有关读取中不匹配数量的证据，因此可以用于筛选出候选区域。

为了有效地使用多路合并，我们的算法依赖于另一种思想，即连续得分约束。对于给定的读取，我们的程序旨在报告“最佳”对齐，即基于不匹配的得分最低的对齐，加上indel或拼接的开口间隙惩罚。因此，我们的搜索过程受到不断增加的分数水平的限制K（K），从开始K（K）=0表示完全匹配，并以成功对齐结束K（K）或用户指定的最高得分水平。除了找到最佳对齐之外，受限搜索过程还可以通过在第一个或最佳的生成对齐的分数级别之外的连续分数级别上继续搜索来找到次优对齐。我们的算法还可以通过搜索给定的分数水平并报告所有结果，找到一组到给定分数水平的详尽对齐。

取决于分数限制K（K）和读取长度我多途径合并过程可以用两种不同的方法来生成和过滤基因组区域。对于低值K（K）不涉及或涉及与我，我们根据生成集寡聚物，根据计数属于q个-支持该地区的mer。对于更高级别的K（K）涉及更多不匹配，我们应用基于成套寡聚物，根据图案属于q个-支持该地区的mer。基于计数和模式的标准都为读取或部分读取中出现的不匹配数量提供了下限。如果下限超过给定的分数限制K（K）在允许的失配中，读取可能会被过滤掉，不需要根据基因组进行验证，以确定失配的实际数量。

哈希表相对较大，如果对每个重叠的寡聚体进行索引，则需要12 GB来表示一个人类大小的基因组。因此，SOAP（Li等。,2008年b)处理人类大小的基因组需要14GB的内存。虽然现代计算机通常有足够的物理内存来查询如此大的哈希表，但较小的数据结构可以通过更有效地使用内存分页和缓存资源来加速程序。我们可以通过对表中索引的基因组低聚物进行采样来减小散列表的大小。在我们的程序中，我们在基因组中每3nt索引12mers，这将人类基因组散列表的大小减少到4GB。因此，我们的算法被设计为使用每3 nt采样一次的哈希表，并且仍然能够像索引每个重叠的低聚物一样实现完全灵敏度。

哈希表索引方案可以扩展到在SNP耐受比对中同样好地比对主要和次要等位基因。（为了便于讨论，我们将参考序列中的等位基因称为“主要”，将其相应的替代版本称为“次要”，而不管它们在群体中的实际频率如何。）因为哈希表代表了q个-mer片段，它可以以一种相对简单的方式表示所有主要和次要等位基因组合的巨大空间。

为了构建一个SNP耐受性哈希表，我们扫描基因组并处理每个采样的基因组q个-包含一个或多个SNP的单体，通过产生包含在其中的每个可能的主要和次要等位基因组合，并为每个产生的基因位置复制该基因组位置q个-梅尔。最后，我们为每个人重新排序位置列表q个-梅尔(图2B） ●●●●。在这个散列表中查找q个-给定基因组位置的mer都将包含所需位置。我们的经验表明，允许SNP的哈希表只比原始哈希表稍大。当我们将来自dbSNP版本129的1200万个SNP合并到人类基因组版本36中时，哈希表的大小从3.8 GB增加到4.0 GB。我们的构造算法要求计算机有足够的内存来存储哈希表，因此人类大小的基因组需要4GB的内存。中讨论了以SNP容忍的方式进行验证第2.4节.

2.2跨度集生成和过滤

跨越集是覆盖读取的最小12毫秒集(图3). 这种结构利用了鸽子洞原理，即不支持12毫秒的数量，即在相应的位置列表中没有包含给定位置的12毫秒，为读取中的不匹配数量提供了一个下限。然而，由于我们在散列表中使用采样，这一鸽子洞原理的实现变得复杂，这就产生了相对于采样的基因组12mers的对齐读取阶段的不确定性。因此，程序必须构造六个生成集，一个生成集用于0、1或2 nt在正向和反向补码方向的移位(图3A） ●●●●。此外，采样哈希表导致基因组位置信息仅以3 nt的间隔可用，从而导致超过读取边界的12毫秒的信息不完整。为了处理这种情况，程序通过替换悬垂位置中所有可能的核苷酸，并取产生的位置列表的并集，来计算悬垂在读取结束时1或2 nt的12个mer的位置列表(图3B） ●●●●。

另一个复杂的问题是，跨越集通常包含12个重叠的节点。为了解决这个问题，我们认为一对重叠的12mer是生成集中的单个“元素”，位置列表等于两个组成位置列表的交集(图3C） ●●●●。由此产生的元素集是不重叠的，因此鸽子洞原理现在适用于k个非支持元素意味着k个不匹配，如果k个>K（K）。对于重叠的一对12mm，可以有多个选择来创建单个元素；我们的程序启发式地选择位置列表最长的12mer或位置列表的并集作为重叠的位置，因为该12mer上的交集操作可能会从后续考虑中消除最大数量的位置。

虽然我们可以使用所有生成集元素来生成候选，然后继续验证步骤，但如果我们指定一些元素来生成候选者，我们可以使算法更快(图3D） 并为单独的过滤步骤保留其他步骤(图3E） ●●●●。这种分工旨在减少O（运行）(我日志n个)基于堆的优先级队列的复杂性，在我。如果我们检查长度的排序列表我_我对于过滤步骤中给定位置的存在，这可以在对数时间内完成O（运行）（日志我_我)通过二进制搜索过程。因此，我们的方法对一些位置列表（生成元素）执行基于堆的合并，并统计支持每个结果候选区域的元素数。如果此计数足够高K（K）或者不支持的总元素更少，则候选区域将执行筛选步骤，检查每个筛选元素是否支持。如果超过K（K）总元素显示不支持；否则，该区域将执行验证步骤，以确定不匹配的实际数量。

我们以各种方式高效地实现了生成集方法。首先，程序选择位置列表最短的元素作为生成元素，选择最长的元素作为过滤元素，因为虽然生成元素中的每个位置都必须处理，但只需要处理过滤元素中的一些位置。其次，我们在每个过滤元素上维护一个指针，只有当我们使用快速二进制搜索（Hwang和Lin，1980). 第三，涉及位置列表的并集或交集的过滤元素不需要显式计算这些集合操作，而是可以用它们的组成位置列表来表示，并在需要时通过执行适当的析取或合取搜索来支持检查。

分配N个生成和筛选目的之间的总元素取决于约束分数级别K（K）允许的不匹配。至少(K（K）+1）当K（K）其他发电元件则不然。我们从经验上发现K（K）>1、分配效率更高(K（K）+2）用于生成目的的元素，因为对两个支持元素的要求大大减少了生成的候选区域的数量。

因为生成集方法至少需要(K（K）+2）生成元件[或(K（K）+1）对于精确和一个不匹配约束]，它只能用于检测相对于读取长度的有限数量的不匹配我，它限制了元素的总数N个。当我=10（mod 12），因此N个≤⌊(我+ 2)/12⌋. 因此(K（K）+ 1)<N个或(K（K）+ 2)<N个指示生成集方法可以应用于约束级别K（K）什么时候K（K）=14≤时为0我≤21;K（K）22≤1≤我≤33; 和K（K）≤⌊(我+2） /12⌋-2用于我≥34.

2.3成套生成和过滤

为了处理比生成集方法检测到的更多的失配，我们采用了一种基于重叠12mm的完整集的策略。此完整集方法适用于任何约束级别K（K）在允许的不匹配中，只要读取区域和候选区域有14个连续匹配（12个月的相位偏差多达2个nt）。连续14场比赛的一个充分条件是失配次数≤⌊我/14⌋−1. 在这种程度的失配情况下，GSNAP是一种详尽的算法，这意味着它可以确保识别和报告基因组中存在如此多失配的所有可用比对。

候选项是通过在单个正向补码和单个反向补码过程中对所有读取位置的位置列表进行多路合并而生成的，并跟踪支持每个候选区域的12毫秒的读取位置。支持12毫秒的模式为读取中的不匹配提供了一个下限。如果支撑12mm的读数位置由Δ隔开第页，则它们之间的最小失配数为⌊（Δ第页+6)/12⌋ (图4A） ●●●●。在整个读取过程中，我们可以在基于模式的下限计算中求和这些下限(图4B） ●●●●。具体来说，如果读取长度我具有在读取位置支持12毫米的模式第页_我,我=1,…,n个，失配的下限是∑_我=0^n个⌊(第页_我+1−第页_我+6） /12⌋，其中第页₀=−3和第页_n个+1=我−9.

为了使全集方法更有效，我们注意到合并过程必须处理每个位置列表中的每个位置，因此，由于位置列表非常长，无法帮助定位读取，因此可以将速度减慢12毫秒。我们可以通过忽略这些非特定的12毫秒来提高效率，目前定义为位置列表长度大于平均位置列表长度的10倍。必须相应地修改下限公式，以补偿缺失的12毫米，基本上是假设它们是支持的。如果非特异性或重复的核苷酸模式是绘制读取图所必需的，则此策略可能无法对齐读取或部分读取。为了成功地对齐这些读数，程序提供了一个贪婪策略选项，在该策略中，最初忽略非特定或重复的12毫秒，如果找不到对齐，则随后包括在内。

2.4候选区域验证

生成的候选区域在过滤过程中幸存下来，它们的不匹配数量有一个确定的下限。确定不匹配的实际数量，并验证其不超过分数限制K（K），我们通过将读数与区域对齐来检查这些区域。为了减少内存需求，我们以压缩格式存储基因组，该格式是在预索引过程中除了哈希表之外创建的(图2C） ，并针对该基因组的压缩版本进行验证。为了索引参考序列，压缩的基因组包含主要等位基因。为了索引参考空间，GSNAP还访问包含次要等位基因的第二个压缩基因组。我们在GMAP论文中描述的压缩基因组格式将每个核苷酸存储在3位中。基因组的每个32nt块由三个32位单词表示。前两个单词分别使用2位表示核苷酸，而剩下的单词使用32位作为标记。对于主要等位基因，标记指示基因组位置是否有未知或不明确的核苷酸，不能用a、C、G或T表示。对于次要等位基因基因组，标记指示基因位置是否有SNP。

在位级别执行验证。该程序不是解压缩基因组，而是压缩读数并将其移位，以匹配候选区域的基因组坐标。然后使用异或函数按位组合压缩的读取和基因组，并减少相邻的位对以产生包含失配位置的位向量。对于参照空间的比对，读取与小等位基因基因组进行类似的位组合，并且使用逻辑和函数组合两个不匹配结果。因此，只有当读取的等位基因与主要和次要等位基因都不同时，SNP才会发生错配。可以进一步分析不匹配结果，以统计不匹配的总数，或从读取的左侧或右侧报告其位置。GSNAP将使用内置的按位函数人口统计,里兹（计数前导零）和ctz公司如果这些任务在给定的机器上可用，则使用它们（计数尾随的零）；如果它们不可用，则将使用其自己的等效位函数。然而，我们的测试表明，内置函数的速度只提高了1-2%，部分原因是生成和过滤步骤大大限制了必须验证的区域数量。

2.5检测插入和删除

我们的程序可以检测包含单个插入或删除的对齐，不匹配的最大值可达用户指定的最大值。相对于无间隙对齐，可以使用用户定义的惩罚对索引对齐进行惩罚G公司因此，随着程序施加越来越强的约束级别K（K）对于允许的不匹配，它还通过施加约束级别来搜索indel对齐(K（K）−G公司)允许与indel不匹配。

索引是使用两种算法检测的，一种是在短读取的中间检测索引（在第一个和最后一个14mer之间），另一种是检测末尾的索引（在第一个或最后一个14 mer内）。末端索引被限制为具有无错配且足够长（由用户指定）的远端短片段，以可靠地确定其对齐。对于这两种方法，执行完整集方法的合并步骤，但不执行过滤计算，以生成与读取内容具有14个或更多连续匹配的所有候选区域。

检测中间索引的方法寻找一对候选区域，这些区域在允许的最大删除和插入大小内共同定位(图5A） ●●●●。基于位置的下限计算可以应用于候选对的两端，以筛选出不匹配过多的对。在核苷酸水平上对其余候选对进行验证，以确定对中每个成员中不匹配的位置，然后可以对其进行分析，以确定约束水平内是否存在可能的交叉区域K（K）允许的不匹配(图5B） ●●●●。中间indel算法即使对于长indel也是有效的，因为两端之间的基因组距离指定了间隙大小，并允许有效验证不匹配，而无需诉诸动态编程算法。

端点索引的检测也依赖于完整集方法生成的候选区域，但会过滤单个候选区域而不是对。使用基于位置的下限计算的变体筛选候选项，该计算忽略读取的第一个或最后一个14mer，该14mer由结束索引设置为不支持。通过此筛选步骤的候选基因将根据基因组进行验证，以计算读取的长部分中的不匹配数量。如果不匹配的数量足够低，则在indel的可能末端插入和删除间隙大小范围内测试末端区域。对于每个间隙大小，程序从末尾开始回溯，直到达到第一个不匹配，然后沿着主对角线计算靠近末尾的不匹配总数(图5C） ●●●●。如果远端段足够长，并且长区域和远端区域的不匹配总和足够低，则检测到末端indel。

2.6检测拼接接头

GSNAP可以比对涉及已知或新剪接位点的跨外显子-外显子连接的转录读数。对于已知的剪接位点，该程序依赖于用户提供的一组剪接位点。剪接位点属于四类之一：正基因组链上的供体和受体，负基因组链上供体和接受体。新剪接位点的识别由概率模型辅助，该模型目前采用最大熵模型（Yeo和Burge，2004)它使用剪接位点附近的核苷酸频率来区分真剪接位点和假剪接位点。

我们使用两种方法检测拼接接头，一种用于短距离拼接，另一种用于长距离拼接。短距离剪接涉及位于同一染色体上的两个剪接位点，受体位点位于供体位点的下游，在用户指定的参数内（默认值为200 000 nt）。除了候选区域之间允许的距离要长得多之外，可以使用与中间删除类似的方法检测短距离拼接接头(图5B） ●●●●。与中间索引检测一样，两个区域中的失配位置决定了交叉区域是否存在允许数量的失配(K（K）−S公司)，其中S公司是接头的开口间隙惩罚。在这个交叉区域中搜索已知或有剪接位点模型支持的供体和受体剪接位点，概率足够高。所需的概率分数取决于外显子区域中可用于比对的短阅读序列的长度。当对齐的外显子序列较短，约为12-20 nt时，需要相对较高的概率分数。但当对齐的外显子序列足够长，超过35 nt时，只需要内含子末端的预期二核苷酸。

对于长距离剪接，概率分数也用于帮助找到新的剪接位点，尽管对于给定长度的对齐序列，所需的概率分数更高，以补偿整个基因组上更大的搜索空间。为了检测长距离拼接的情况，该程序在由约束级别限定的区域内识别单个候选区域内的已知或新的拼接端K（K）允许的不匹配(图5D） 。如果具有供体剪接位点和受体剪接位点的候选区域在读取时具有相同的断点，并且具有可接受的总失配数，则对其进行配对。

主要位于剪接接头一端的读数可能在远端的序列太少，无法识别其他外显子。如果一端有足够的序列来确定一个拼接位点，而另一端没有足够的序列用于另一个位点，我们的程序仍然可以将这种比对报告为部分拼接或“半内含子”比对。

2.7对齐成对读取

GSNAP可以对齐配对读码，这些读码是在核苷酸片段两端测序时产生的。成对读取也可以通过使用短链接器循环10000个或更多碱基的长片段来生成，然后在链接器外部进行剪切，以在长片段的两端进行读取。我们的算法试图找到一对和谐的最佳映射变体，这意味着它们在用户定义的预期基因组距离范围内，并且它们的链方向是一致的。因此，GSNAP将支持涉及一端或两端次优对齐的协调解决方案，即使可以为每一端单独找到更好的对齐。

为了将成对的路线放在一起考虑，程序会连续施加更强的约束级别，并尝试在给定的约束级别上对齐每一端。对于成对对齐，读取的两端都有助于获得总分，因此在给定的约束级别上K（K），程序必须将每个端点的对齐累积到该水平。在每个约束级别上，程序会尝试配对到目前为止在每一端找到的一组详尽的对齐，以查看是否有任何配对是一致的。如果是这样的话，最好的一条或多条路线是那些在两端不匹配和惩罚总分最低的路线。如果用户需要次优比对，程序将继续查找超出最佳评分水平的其他比对。如果该算法在没有找到任何一致对的情况下达到了用户定义的最大分数限制，它会报告每一端的最佳单独对齐。

2.8 BS-转化DNA

辅助程序处理现有的参考序列或参考空间散列表以生成两个新的散列表，这两个散列表都表示基因组的正链，其中一个具有C–T替换，另一个具有G–A替换。第二个哈希表容纳从负链读取的数据，其C–T替换在正链上显示为G–A替换。辅助程序将每个被替换的12毫米的位置组合并排序为一个位置列表。

当GSNAP处理一个BS读取时，它在读取过程中对每个12mer执行C–T替换以检查C–T哈希表，并在读取的反向补码中对每个12 mer执行G–a替换以检查G–a哈希表。生成和过滤步骤的行为与之前相同。验证步骤将取代的读数与取代的基因组区域进行比较，以识别失配。在原始读数中对照原始基因组区域进行特殊检查，以确定基因组-T和读数-C之间的不匹配，这些不匹配被C-T取代所掩盖。

3结果

3.1模拟读取

我们将GSNAP与之前研究中已进行基准测试的几个校准程序进行了比较：MAQ版本0.7.1、SOAP版本1.11、Bowtie版本0.9.9.1、BWA版本0.4.9和SOAP2版本2.19。我们生成了36个、70个和100个从人类基因组中统一采样的nt读取（NCBI第36版），并通过随机引入不匹配或索引生成了不同变体类型的数据集。短indels（1–3 nt）要求距离末端至少6 nt，长indels至少14 nt。程序在Linux机器上运行，该机器配有8个双核AMD 8220 Opteron CPU，2.8 GHz和64 GB RAM。所有程序都具有多线程模式，但在我们的测试中是以单线程模式运行的。为了分别研究每个变体，并防止程序搜索次优点击，我们为每个程序提供了足以识别给定变体的参数。

对齐结果(表1)显示程序通常可以在最多三个不匹配的情况下正确对齐读取，SOAP2限制为两个不匹配，GSNAP和SOAP的36 nt读取和SOAP的70 nt读取中观察到的未命中。GSNAP在36 nt读取中两个和三个失配的缺失是指失配数量超过⌊保证条件的情况我/14⌋−1表示穷尽性，主要是因为错配的间隔足够均匀，以防止连续延伸14 nt以匹配读取和基因组。Bowtie和SOAP2最快确定了精确比对和单错配比对。BWA和GSNAP以可比较的速度识别70和100 nt读取中两个和三个不匹配的对齐，而GSNAP在识别四个或五个不匹配时速度最快。对于短indels，GSNAP显示出完美的灵敏度，而BWA和SOAP显示出高达5%的漏检率。SOAP2无法在我们的数据集中检测到单端读取中的索引，尽管它可以在成对读取中检测到短索引。在除一个数据集外的所有数据集中，GSNAP在识别短indels方面速度最快。对于长索引，GSNAP是唯一能够检测到对齐的程序，但SOAP也可以以较慢的速度检测36 nt读取中的长索引。GSNAP检测长indels的速度与短indels相当，并且显示出0.1%的漏检率，这都是由于indels一端的重复区域造成的，程序设计为忽略该区域。

我们使用Valgrind Massif工具（Nethercote和Steward，2007). GSNAP在完全匹配的数据集上使用了86 MB的峰值，在不匹配的数据集上使用了101 MB，在indel数据集上使用了170 MB，内存使用情况因读取而异。然而，这些值并不测量GSNAP对其大小为3.8 GB的哈希表位置文件和大小为1.15 GB的压缩基因组使用的内存，这些数据在主机上可用时使用内存映射进行访问。在内存映射中，当有足够的物理内存来保存索引文件的相关部分时，我们的程序将以最快的速度运行。因此，为了在人类大小的基因组上实现最佳性能，GSNAP应该可以访问5 GB的物理内存，尽管如果可用内存较少，程序仍然可以运行，尽管速度会更慢。内存映射的一个优点是，GSNAP的多个实例可以同时在同一台计算机上运行，并共享映射到索引文件的系统内存，而无需每个进程分别分配该内存。为了进行比较，BWA在所有数据集上使用了2.2 GB；Bowtie在精确匹配的数据集上使用1.1 GB，在不匹配的数据集中使用2.2 GB；MAQ在所有数据集上使用了302 MB；SOAP在所有数据集上使用了14GB。无法确定SOAP2的内存使用情况，因为源代码不可用，并且编译二进制程序时没有内存分析所需的标志。

我们还评估了GSNAP在模拟阅读中检测基因内和基因间外显子-外显子连接的能力，使用RefSeq中的已知剪接位点或仅使用新的剪接位点检测（这些测试是在使用绝对概率阈值的早期版本的GSNAP上进行的，而不是在当前版本中使用的滑动尺度）。模拟读取基于RefSeq剪接位点，GSNAP在获得该信息时达到了完美的灵敏度，但在仅依赖概率剪接位点模型时，遗漏了0.1%的基因内事件和5%的基因间事件。缺失剪接事件是由于已知剪接位点的模型得分低于默认概率阈值0.90。基因内和基因间剪接检测的敏感性差异是由于GSNAP中短距离剪接的标准不同，后者对第二个剪接位点使用较低的概率阈值（默认值为0.50）。对于36、70和100 nt的基因内读取，10万次读取的数据集的运行时间分别为48、114和183秒，对于基因间读取，运行时间分别是122、211和287秒。这些运行时间比不匹配和索引的基准测试更快，因为拼接读取是从基因组的编码部分生成的，重复性较小。

3.2转录阅读

我们测量了GSNAP对复杂变异检测和SNP耐受性的影响，这些影响来自MAQC（微阵列质量控制）项目（Canales）中使用的通用人类参考RNA（UHR，Stratagene目录号740000）的实际数据等。,2006)并由Illumina在Solexa基因组分析仪上进行分析。从这个数据集中，我们在50 nt个读取的数据集中均匀采样了10万个读取。与模拟读取不同，实际读取缺少关于其原始基因组位置的信息，因此我们使用对齐率来确定程序的性能，对齐率是每个程序在各种不匹配、剪接或索引设置下可以对齐的读取的百分比。我们针对NCBI人类基因组版本36测试了GSNAP，还针对覆盖了来自dbSNP版本129的1200万SNP的参考空间进行了测试。

两个失配的对准率为70%，三个失配为74%(表2). 当允许程序识别更复杂的变体时，校准产量增加。在GSNAP中添加包含已知剪接位点的剪接增加了8-9%的对齐产量，而进一步添加新的剪接位点又增加了0.3-0.6%的对齐产量。允许使用indels将对齐收益率再提高1%。

引入SNP公差仅导致对齐产量略有增加。然而，7-8%的比对结果在某种程度上受到SNP的影响。SNP耐受性与之前只有0.4-0.5%的病例未发现SNP耐受的情况一致。在5%的病例中，已知的SNPs在原始位置之外显示了给定读数的额外基因组位置，从而产生了原始结果的超集。在0.2-0.4%的病例中，SNP耐受性产生了基因组位置的子集，这意味着原始比对中的一些错配可以通过已知SNP解决。在1-2%的病例中，原始比对中的所有错配均位于已知的SNP位置，使基因组位置保持不变，但允许核苷酸差异被解释为与次要等位基因的匹配，而不是错配。在一小部分病例中，SNP耐受性给出了与原始结果显著不同的结果。

3.3 BS-转换读取

我们使用了来自第3.1节具有0、1和2nt的错配，并且以95%的概率用胸腺嘧啶取代胞嘧啶，忽略序列上下文，例如真核生物中的非岛CG二核苷酸（Goll和Bestor，2005)或植物中的CG、CHG和CHH模式（Cao和Jacobsen，2002)其中甲基胞嘧啶出现频率更高。我们将原始读取结果与标准版本的GSNAP进行比对，并将甲基化标记打开的替换读取结果进行比对。结果表明，胸腺嘧啶替代对GSNAP识别原始基因组位置的能力影响很小，因为比对某些读取结果的模糊性增加(表3). 在完全匹配的数据集中，在基因组中给出非唯一位置的额外读取的比例为36 nt读取的5.5%、70 nt读取的2.6%和100 nt读取的1.3%，在有一个或两个不匹配的数据集中，比例稍高。

4讨论

本文描述的方法扩大了可在读取中检测到的变体的范围，因此应提高下一代测序数据的效用。识别广泛变异的能力还应通过识别变异读取的正确基因组起源来提高绘图准确性。同样，我们程序中实现的SNP耐受特性应该有助于解决某些基因组区域的映射。此功能的实用性不仅应通过7-8%的受影响读取来衡量，还应通过其在后续分析管道中对做出正确生物推断的贡献来衡量。其他研究人员（Manske和Kwiatkowski，2009)也发现了SNP耐受性比对有助于阅读与次要等位基因的比对的情况。

我们开发了一种算法，以满足短读序列分析在检测复杂变体和剪接时的速度和灵敏度方面的特定需求。我们算法的优势在于它的连续约束搜索策略，通过合并整个读取过程中寡聚体的位置列表，并使用基于计数或模式的下限计算对其进行过滤，从而生成候选基因组区域。我们的搜索过程在低聚物水平上运行，这与基于BWT的程序中用于识别不匹配和短索引的核苷酸水平回溯过程不同。

通过筛选候选区域集，我们的交集过程比BLAT和其他基于种子的比对程序中使用的种子和扩展策略的效率有了显著提高，这些程序基于单个区域找到基因组区域q个-mer，然后在耗时的验证步骤中测试每个区域。一些程序，如ELAND，进一步将种子限制在短读的开头。种子的使用可以是一种高效的启发式方法，而BWA可以通过使用种子模式来运行得更快，这种模式允许在简短读取的初始部分出现一定数量的不匹配。我们的方法可以被认为是在整个短阅读过程中同时尝试所有可能的种子，因此具有不偏爱阅读的一部分而不偏爱另一部分的特点。

我们的交叉程序也为q个-克程序，用于SHRiMP（隆隆声等。,2009)和RazerS，它使用滑动窗口扫描整个基因组并计数q个-mers在箱子中寻找候选基因组区域（拉斯穆森等。,2006). 这个q个-gram过程允许一次读取两个或多个索引，这是GSNAP当前不允许的，尽管我们的算法可以修改以识别它们。另一个区别是q个-gram程序将比对差异定义为编辑距离，其中间隙中的每个核苷酸都算作差异，所以indel越长，距离越远。相比之下，GSNAP在得分indel对齐时仅使用了一个开口间隙惩罚，因此它可以更容易地识别长indel。

我们的计划也不同于几个校准计划，包括MAQ、RMAP、SHRiMP、Bowtie、BWA和SOAP2，它们能够使用质量分数对校准进行排名。尽管质量分数可以应用于我们算法的验证步骤，但我们仍不清楚如何在质量分数和对齐结果之间进行最佳权衡，例如，如何在一个高质量分数的一个不匹配对齐和一个低质量分数的两个不匹配的对齐之间进行选择。ABI SOLiD技术产生的颜色空间读取需要对我们的算法进行一些扩展，我们正在努力在我们的程序中实现这一功能。

尽管我们的结果和经验表明，我们的程序对于分析下一代测序数据具有实用价值，但我们的研究仍在继续。特别是，更长的读数将需要更大的比对灵活性，并可能需要加强更通用的cDNA-基因组比对程序，如我们的GMAP程序。未来的生物研究应受益于具有多样性的生物信息学方法和程序，以满足序列分析的各种需求。

致谢

我们感谢同事Colin Watanabe的宝贵讨论，感谢Sekar Seshagiri在下一代测序项目上的合作。我们感谢国家基因组资源中心的Andrew Farmer和Ernie Retzel对我们项目早期版本的反馈。我们还感谢Illumina的Irina Khrebtukova和Gary Schroth访问转录读取数据。

利益冲突：未声明。

参考文献

作者注释