STAR:超快通用RNA-seq对准器

摘要

动机:高通量RNA-seq数据的精确比对是一个具有挑战性但尚未解决的问题，因为测序技术的不连续转录本结构、相对较短的读取长度和不断增加的通量。目前可用的RNA-seq校准器存在高映射错误率、低映射速度、读取长度限制和映射偏差。

结果:为了对齐我们的大型(> 800亿读取)ENCODE转录组RNA-seq数据集，我们开发了基于先前描述的RNA-seq对齐算法的拼接转录本比对参考(STAR)软件，该算法使用在未压缩后缀数组中顺序最大可映射种子搜索，然后进行种子聚类和拼接过程。STAR在定位速度上比其他比对器高出>50倍，在一个中等的12核服务器上，每小时对人类基因组进行5.5亿次2 × 76 bp配对端读取，同时提高了比对灵敏度和精度。除了不偏不倚新创STAR可以发现非标准剪接和嵌合(融合)转录本，也能够映射全长RNA序列。我们利用Roche 454对逆转录聚合酶链反应扩增子进行测序，实验验证了1960个新的基因间剪接连接，成功率为80-90%，证实了STAR作图策略的高精度。

可用性和实现:STAR是作为一个独立的c++代码实现的。STAR是在GPLv3许可下发布的免费开源软件，可以从http://code.google.com/p/rna-star/．

联系人:dobin@cshl.edu．

1介绍

虽然基因组是由线性有序的核酸序列组成，但真核细胞通常通过拼接不连续的外显子来重组转录组中的信息，以创建成熟的转录本(Hastings和Krainer, 2001)．这些剪接rna的检测和表征一直是正常和疾病细胞状态下基因组功能分析的关键焦点。测序技术的最新进展使得单核苷酸水平的转录组分析几乎成为常规。然而，由这种高通量测序实验产生的数以亿计的短(36 nt)到中(200 nt)长度序列(reads)对拼接转录本的检测和表征提出了独特的挑战。有两个关键任务使得这些分析需要大量的计算。第一个任务是对包含由基因组变异和测序错误引起的不匹配、插入和删除的reads进行精确校准。第二项任务涉及映射来自非连续基因组区域的序列，这些区域包含拼接序列模块，这些拼接序列模块连接在一起形成拼接rna。虽然第一个任务是与DNA重测序工作共享的，但第二个任务对RNA-seq是特定的和至关重要的，因为它提供了重建全部剪接RNA分子所需的连通性信息。由于存在多个相同或相关的基因组序列副本，这些序列本身也被转录，使得精确的定位变得困难，从而进一步加剧了这些校准挑战。

各种序列比对算法最近被开发来应对这些挑战(非盟et al。, 2010年；德善et al。, 2008年；格兰特et al。, 2011年；汉et al。, 2011年；杰尔et al。, 2009年；王et al。, 2010年；Wu和Nacu, 2010；张et al。, 2012年)．然而，这些算法的应用需要在映射精度(灵敏度和精度)和计算资源(运行时和磁盘空间)方面做出妥协(格兰特et al。, 2011年)．随着当前测序技术的进步，计算组件日益成为吞吐量的瓶颈。高映射速度对于大型财团的工作尤其重要，例如ENCODE (http://www.genome.gov/encode/)，不断生成大量的测序数据。

此外，大多数引用的算法设计用于处理相对较短的读取(通常≤200个碱基)，不适用于对齐新兴的第三代测序技术生成的长读取序列(Flusberget al。, 2010年；Rothberget al。, 2011年)．较长的读取序列，理想情况下达到RNA分子的完整长度，通过提供更完整的RNA连通性信息，在增强转录组研究方面具有巨大潜力。

本报告描述了一种名为“拼接转录本对参考(STAR)”的对齐算法，该算法旨在专门解决RNA-seq数据映射的许多挑战，并使用了一种新的拼接对齐策略。我们进行了高通量验证实验，证实了STAR在检测新型剪接连接方面的精度。STAR的高绘制速度和准确性对于分析大型ENCODE转录组至关重要(Djebaliet al。, 2012年)数据集(> 800亿Illumina读取)。我们还证明了STAR在精确对齐第三代测序技术中出现的长序列(几千碱基)读取方面的潜力。

2算法

许多先前描述的RNA-seq比对器被开发为连续(DNA)短读映射器的扩展，用于将短读对齐到剪接连接数据库或将分裂读部分连续对齐到参考基因组，或两者的组合。与这些方法相比，STAR被设计成将非连续序列直接对准参考基因组。STAR算法包括两个主要步骤:种子搜索步骤和聚类/拼接/评分步骤。

2.1种子搜索

寻找STAR种子阶段的中心思想是对最大可映射前缀(MMP的)．MMP的类似于大规模基因组比对工具使用的最大精确(唯一)匹配概念Delcheret al。, 1999年，2002；库尔茨et al。)和淡紫色(亲爱的et al。, 2004年，2010)。给定一个读序列R，读取位置我和一个参考基因组序列G,MMP(右，我，G)定义为最长子字符串(R_我，R_我+1，…，R_我+MML−1的一个或多个子字符串G,在那里MML是最大可映射长度。我们将用一个简单的例子来解释这个概念，这个例子是一个包含单个拼接结且没有错配的读操作(图1a).在第一步中，算法找到MMP的从读取的第一个碱基开始。由于本例中的reads包含一个剪接结，因此它不能连续地映射到基因组，因此第一粒种子将被映射到一个供体剪接位点。接下来,MMP的重复搜索读取的未映射部分，在这种情况下，将映射到一个受体拼接位点。注意，这种顺序的应用MMP的只搜索读取的未映射部分使得STAR算法非常快，并将其与Mummer和MAUVE区别开来，后者可以找到所有可能的最大精确匹配。这种方法代表了一种在读取序列中找到剪接结点精确位置的自然方法，并且比在分割-读取方法中使用的读取序列的任意分割更有优势。拼接连接是在单个对齐通道中检测到的，没有任何连接先天的了解拼接节点的位置或属性，并且没有连接数据库所需的初步连续对齐通道。的MMP的在STAR中通过未压缩后缀数组(sa)实现搜索(曼伯和迈尔斯，1993)。值得注意的是,发现MMP的是未压缩sa中标准二进制字符串搜索的固有结果，与全长精确匹配搜索相比，它不需要任何额外的计算工作。SA搜索的二进制性质导致搜索时间与参考基因组长度的有利对数缩放，即使针对大基因组也允许快速搜索。对每个人都有利MMP的SA搜索可以用很少的计算开销找到所有不同的精确基因组匹配，这有助于对映射到多个基因组位点的读取进行精确对齐(“multimapping”读取)。

除了检测剪接结，MMP的在STAR中实现的search可以查找多个不匹配和索引，如图1b.如果MMP的如果存在一个或多个不匹配，则搜索不会到达读取的结束MMP的S将作为可以扩展的基因组中的锚点，允许错配的对准。在某些情况下，扩展过程不能产生良好的基因组比对，这允许识别poly-A尾、文库适配器序列或测序质量差的尾(图1c)。MMP的在读取序列的正向和反向两个方向进行搜索，并且可以从整个读取序列的用户定义搜索起点开始，这便于为接近末端的错误读取找到锚点，并提高了高测序错误率条件下的映射灵敏度。

除了高效MMP的在搜索算法中，未压缩sa也比许多流行的短读对齐器中实现的压缩sa具有显著的速度优势(补充部分1.8)。这种速度优势与未压缩数组增加的内存使用量相权衡，这将在3.3节中进一步评估。

2.2聚类、拼接、评分

在算法的第二阶段，STAR通过将第一阶段中与基因组对齐的所有种子拼接在一起，构建整个read序列的对齐。首先，这些种子通过靠近一组选定的“锚定”种子而聚集在一起。我们发现，选择锚点的最佳程序是通过限制锚点对准的基因组位点的数量。所有在用户定义的基因组窗口内围绕锚点映射的种子被缝合在一起，假设一个局部线性转录模型。基因组窗口的大小决定了剪接序列的最大内含子大小。一种节约的动态规划算法(见补充部分1.5(详细信息)用于缝制每对种子，允许任意数量的不匹配，但只有一次插入或删除(间隙)。

重要的是，配对端RNA-seq读对的配偶的种子是同时聚集和缝合的，每个配对端读对表示为单个序列，允许配偶的内端之间可能存在基因组间隙或重叠。这是使用配对端信息的一种有原则的方式，因为它更好地反映了配对端读取的性质，即配对是相同序列的片段(末端)。这种方法增加了算法的灵敏度，因为只有来自一个同伴的一个正确的锚足以准确地对齐整个读取。

如果一个基因组窗口内的比对没有覆盖整个读取序列，STAR将尝试找到两个或更多覆盖整个读取序列的窗口，从而导致嵌合比对，读取的不同部分映射到远端基因组位点，或不同的染色体，或不同的链(补充图S1)．STAR可以发现嵌合排列，其中配偶之间是嵌合的，其中一个嵌合连接位于两个配偶之间RNA分子的未测序部分。STAR还可以发现其中一个或两个配偶内部嵌合对齐的嵌合对齐，从而精确定位基因组中嵌合连接的精确位置。从K562红白血病细胞系中给出了一个BCR-ABL融合转录本检测的例子补充部分1.7 (补充图S2)．

拼接由局部对齐评分方案指导，对匹配、不匹配、插入、删除和拼接连接间隙使用用户定义的分数(惩罚)，允许对对齐质量和等级进行定量评估(参见补充部分详见1.4)。得分最高的缝合组合被选为读取的最佳对齐。对于multimapping读取，将报告分数低于最高分数的某个用户定义范围内的所有对齐。

尽管是连续的MMP的搜索只找到与基因组完全匹配的种子，后续的拼接程序能够将大量不匹配的reads、indels和splice连接进行对齐，并随着读取长度的增加而扩展。随着第三代测序技术(如Pacific Biosciences或Ion Torrent)的出现，这一特性变得越来越重要，这些技术可以产生更长的读取时间，出错率也更高。

3的结果

3.1模拟RNA-seq数据的性能

首先，我们使用模拟数据来评估STAR的性能，并将其与其他RNA-seq绘制器进行比较。模拟允许精确计算假阳性和阴性率，尽管用于生成模拟读数的人工误差模型可能不能充分代表实验误差。我们使用了来自最近研究的模拟数据集(格兰特et al。, 2011年)，从小鼠转录组中生成了1000万个2 × 100 nt的Illumina-like读序列，包括注释转录本和人工转录本，错误率较高。引入各种类型的基因组变异和测序错误来模拟真实的RNA-seq数据。

STAR 2.1.3、TopHat2 2.0.0 (杰尔et al。, 2009年)， gsnap 2012-07-03 (Wu和Nacu, 2010)、朗姆酒1.11 (格兰特et al。, 2011年)和MapSplice 1.15.2 (王et al。, 2010年)在标记为“SIM1-TEST2”的模拟数据集(格兰特et al。, 2011)。因为TopHat2 2.0.0版本代表了TopHat校准器的主要新发展，它还没有经过同行评审，我们还与之前的TopHat版本1.4进行了比较。我们发现新版本产生了更好的精度和更快的映射速度(补充部分2.1和补充数据)．所有的对准器都在新创模式，即不使用基因/转录本注释。最大不匹配数设置为每对末端读取10个，最小/最大内含子大小设置为20 b/500 kb (补充部分2，以获取更多信息)。请注意，在映射器与默认参数之间运行比较是一种合理且普遍接受的做法，因为所有考虑的校准器默认情况下都针对哺乳动物基因组和最近的RNA-seq数据进行了优化。

将得到的比对结果与模拟reads的真实基因组起源进行比较，并使用由格兰特et al。(2011)。ROC曲线(图2)通过每个结点映射的读取数给出的检测(辨别)阈值来计算，即对于每个校准器，在ROC曲线上的每个点只选择至少有N个读取支持的节点，N从1(最低阈值)到100(高阈值)不等。所有对准器在检测阈值较高时均表现出理想的陡峭ROC曲线。在每个结1个读取的最低检测阈值下，STAR表现出最低的假阳性率，同时实现高灵敏度。补充图S5显示了对低错误率模拟数据集的相同分析，得出了类似的结论。

3.2实验RNA-seq数据的性能

为了评估RNA-seq绘图器在实验RNA-seq数据上的性能，运行STAR、TopHat2、GSNAP、RUM和MapSplice(参见补充部分在ENCODE长RNA-seq数据集(K562全细胞A +样本，1 Illumina GAIIx车道4000万2 × 76读取)。STAR和GSNAP对齐的读取百分比最大(均为94%)，其次是RUM (86%)， MapSplice(85%)和TopHat2(71%)。

与Gencode 7相比，不同的拼接结检测精度指标(哈罗公学et al。, 2012年)注释被绘制在图3a - c作为检测阈值的函数，定义为每个结的RNA-seq读取的最小数量。尽管所有对齐器都检测到相同数量的带注释的连接(图3A，实线)，在检测到的未注释连接的数量(图3A，虚线)。在所有检测到的连接中绘制未注释的连接的百分比图3B作为检测阈值的函数。因为所有对齐器对带注释的连接显示出相似的敏感性，所以所有检测到的带注释的连接的比例可以作为精度的替代品。STAR, RUM和TopHat2的表现相似，而GSNAP在较低的检测阈值下表现出较低的精度，MapSplice表现出不寻常的非单调和非饱和行为，这在中也被注意到张et al。(2012)。伪roc曲线，即检测到的标注节点的比例(伪敏感性)与检测到的未标注节点的比例(伪假阳性率)图3c.所有校准器(MapSplice除外)在检测阈值较高时的性能相似。

因为许多未注释的连接代表了真正的新剪接事件，而不是假阳性，所有检测到的连接中未注释的百分比并不能准确地代表假阳性率。为了更准确地估计假阳性率，我们采用了另一种常用的方法(张et al。, 2012)并绘制(图3D)至少两个映射器检测到的连接数(伪真阳性)和每个映射器单独检测到的连接数(伪假阳性)。STAR比对产生最低的伪假阳性率，即唯一检测到的连接比例最低(图3E)，同时实现了类中第二的伪灵敏度(图3f). GSNAP以较高的伪假阳性率为代价，表现出最高的伪敏感性。这些结果在定性上与对准器在模拟数据上的表现一致，而定量上的差异可能归因于真实误差和模拟误差之间的差异。补充图S6显示了对更短的RNA-seq数据集(2 × 50 b)的相同分析，这表明STAR即使对短读取也保持高灵敏度和精度。

请注意，伪真/假阳性定义是基于这样的假设:只有一个校准器检测到的连接比两个或多个校准器检测到的连接更有可能是假阳性;然而，这些定义并不严谨，因为无法对实验数据进行真假评估。我们还想强调的是，这些比较是针对每种工具的当前版本、默认参数和Illumina测序技术的当前状态进行的。随着测序技术和工具的改进，这些排名可能会改变，必须重新评估。

与其他RNA-seq校准器类似，STAR的默认参数是针对哺乳动物基因组进行优化的。其他物种可能需要对一些对齐参数进行重大修改;特别是，对于内含子较小的生物体，最大和最小内含子大小必须减小。

3.3速度基准

在一台配备了两个6核Intel Xeon cpu X5680@ 3.33 GHz和148 GB RAM(随机访问内存)的服务器上执行速度基准测试。每次运行请求6或12个线程，使用服务器的一半或全部容量。在上一节中描述的~ 4000万2 × 76 Illumina人类RNA-seq数据集上，所有绘图仪都使用默认参数运行。

“墙”时间(即完成映射所需的总运行时间)和RAM使用显示在表1．STAR使用12个线程(服务器的全容量)实现5.5亿2 × 76 Illumina配对端读取每小时的速度，即每处理器每小时4500万配对读取，比第二快的映射器(TopHat2)高出50倍。STAR显示出吞吐量随线程数的线性缩放，当线程数从6增加到12时，每个线程映射速度会损失约10%。

STAR的高绘制速度与RAM的使用是相互权衡的:STAR需要约27 GB的RAM才能对准人类基因组。与所有其他对齐器一样，除了RUM之外，STAR使用的RAM量不会随着线程数量的增加而显著增加，因为SA是在所有线程之间共享的。虽然STAR的RAM要求在几年前非常昂贵，但当第一个短读取校准器被开发出来时，最近半导体技术的进步导致RAM价格大幅下降，现代高性能服务器通常配备RAM >32 GB。STAR可以选择使用稀疏sa，将人类基因组的RAM消耗降低到<16 GB，但在保持校准精度的同时，测绘速度降低了约25%。

3.4实验验证

作为人类转录组特征的一部分，ENCODE (Djebaliet al。, 2012年)， STAR用于绘制从原代人H1ES(胚胎干细胞)和HUVEC(脐静脉内皮细胞)细胞系全细胞提取物中分离的多聚腺苷酸(poly A+)长(>200 nt)转录本。这些rna使用双工特异性核酸酶协议进行测序(Parkhomchuket al。, 2009年)，生成2 × 76 bp的特定于链的读取。

毫不奇怪，未标注的(新)剪接位点比标注的连接显示出更低的丰度水平，这由支持读数的未标注连接数量的显著下降所表明(补充图S7)．由于每个细胞系都是在生物副本中测序的，因此可以根据它们在副本之间的再现性来识别高置信度剪接位点的集合。为了评估检测到的剪接连接的重现性，我们开发了一种非参数不可重现发现率(npIDR)方法，特别适用于RNA-seq数据的离散性质(参见补充材料详细描述)。这种方法类似于ENCODE ChIP-seq实验分析中广泛使用的npIDR概念(Landt编著et al。, 2012年)．补充图S8显示了npIDR = 0.1对每个结的读取计数的依赖性，为选择具有理想再现性水平的读取计数阈值提供了原则性方法。例如，每个连接需要5个交错的读取，以达到0.1的npIDR，即在同一细胞系上具有相同测序深度的另一个实验中再次观察到这些连接的90%的可能性。

实验验证在广泛的RNA-seq读取支持的1920个新剪接连接上进行，均低于或高于npIDR阈值。只有剪接连接映射到Gencode 7基因的基因间或反义位点(哈罗公学et al。, 2012年)被选择进行验证，因为这些连接比注释基因内的连接更有可能是假阳性。高通量验证管道包括目标区域的逆转录聚合酶链反应扩增，然后对汇集产物进行Roche 454测序。逆转录聚合酶链式反应引物设计利用了支持目标连接的配对末端reads的~ 250 nt插入长度，并需要生产300-600 nt长的扩增子。这些扩增子由Roche 454测序仪汇集并测序，以提供长且更有信心的可映射读数，这些读数与BLAT对准基因组。实验协议的详细描述可以在Djebaliet al。(2012)。

我们从H1ES和HUVEC细胞系中选择了1920个基因间和反义剪接连接，包括高(npIDR < 0.1)和低(npIDR > 0.1)可重复的连接。在所有测试的新基因间/反义连接中，至少有5个RNA-seq reads(对应于npIDR < 0.1)支持，约82-89% (H1ES)和84-95% (HUVEC)被至少两个由454 (表2)．值得注意的是，即使只有两个RNA-seq reads支持的候选连接，验证率仍然保持在72% (H1ES)和74% (HUVEC)的高水平。这些结果证实了STAR的拼接检测算法即使对于罕见的新连接也具有很高的精度。

根据验证率(≡VR)，可以估计出错误发现率(FDR)的上限FDR≤1−VR。对于低丰度连接，实验FDR低于由重复之间的差异预测的npIDR:例如，尽管只有两个reads支持的45%的连接是不可重复的(补充图S8)， > 70%的验证成功(表2)．因此，在验证实验不切实际的情况下，npIDR可以作为保守的FDR估计上限。

4讨论

尽管经过几年的不断改进，将不连续的RNA-seq读取序列对准参考基因组还不是一个解决的问题，这既是由于其固有的复杂性，也是由于测序技术的快速转变。已经发现一些严重的问题困扰着以前发表的方法，例如高映射错误率、对齐偏差、对无注释转录本的低敏感性、读取长度的可伸缩性差、每次读取的连接/错配/ indes数量的限制、无法检测非线性转录本(如嵌合rna)，以及关键的低映射吞吐量。

在这项工作中，我们描述了STAR，一种将高通量长和短RNA-seq数据与参考基因组对齐的新算法，旨在克服上述问题。与许多其他RNA-seq绘图仪不同，STAR不是短读DNA绘图仪的扩展，而是作为独立的c++代码开发的。STAR能够在多核系统上运行并行线程，并且随着核数的增加，其工作效率近乎线性扩展。STAR的速度很快:在一台现代化但不太昂贵的12核服务器上，它可以每小时对齐5.5亿个2 × 76 nt的人类基因组，超过所有其他现有的RNA-seq校准器50倍。与此同时，STAR在实验和模拟数据上都表现出比其他RNA-seq校准器更好的校准精度和灵敏度。

所有问题中最主要的一个新创RNA-seq校准器无法准确检测涉及短(< 5-10 nt)序列在连接的供体或受体侧悬垂的剪接事件。这导致剪接事件的显著不足检测，也显著增加了错配率，因为这样的读取可能与一些类似的连续基因组区域不匹配。此外，这种效应也会使排列偏向加工过的假基因，这些假基因在人类基因组中非常丰富。与其他RNA-seq校准器类似，为了缓解这个问题，STAR有一个选项，可以从注释数据库(补充部分4).也可以运行第二个映射通道，为其提供在第一个映射通道中发现的剪接结位点。在这种情况下，STAR将不会发现任何新的连接，而是将拼接的读取与先前检测到的连接上的短悬垂对齐。

为了证明STAR对齐长读取的能力，我们绘制了来自GenBank的长(0.5-5 kb)人类mRNA序列补充部分5)。STAR对齐的准确度与BLAT相近或更高(肯特,2002)一种流行的EST/mRNA校准剂。同时，STAR在比对速度上比BLAT高出两个数量级以上，这对于高通量测序应用非常重要。

算法对长读取的可扩展性表明STAR有潜力作为一种通用比对工具，适用于广泛的新兴测序平台。STAR可以在连续流模式下对读取数据进行对齐，这使得它与新的测序技术兼容，例如牛津纳米孔技术公司最近宣布的一项测序技术。随着测序技术和协议的发展，必须开发新的映射策略，STAR核心算法可以提供一个灵活的框架来解决出现的对齐挑战。

数据访问

GEO: GSE38886 (Roche 454测序)

GEO: GSE30567 (Illumina long RNA-Seq)

资金:本工作由NHGRI (NIH)资助U54HG004557。

利益冲突:未声明。

参考文献

作者指出