从RNA-seq读取到差异表达结果

摘要

许多方法和工具可用于预处理高通量RNA测序数据和检测差异表达。

高通量测序技术目前在生物学中普遍使用。这些技术产生了数以百万计的短序列读取，并定期应用于基因组、表观基因组和转录组。在样本中对稳态RNA进行测序（称为RNA-seq），不受以往技术的许多限制，如微阵列和PCR所需的对生物体先验知识的依赖（见方框1：微阵列和基因表达分析测序的比较）。此外，RNA-seq有望解开转录组中以前难以接近的复杂性，如等位基因特异性表达和新的启动子和亚型[1——4]. 然而，所产生的数据集又大又复杂，解释也不简单。与任何高通量技术一样，分析方法对解释数据至关重要，RNA-seq分析程序也在不断发展。因此，及时回顾当前可用的数据分析方法并评论未来的研究方向。

理解RNA-seq数据取决于感兴趣的科学问题。例如，确定等位基因特异性表达的差异需要准确测定转录单核苷酸多态性（SNP）的流行率[5]. 或者，可以通过在RNA-seq数据中发现新的转录物来检测癌症样本中的融合基因或畸变[6,7]. 在过去的一年中，出现了几种使用RNA-seq数据进行丰度估计的方法[8,9]，检测替代拼接[10——12]，RNA编辑[13]和新的抄本[11,14]. 然而，许多生物学研究的主要目标是样本之间的基因表达谱分析。因此，在这篇综述中，我们重点介绍了检测样本之间基因水平表达差异的可用方法。这种分析特别适用于对照实验，比较同一组织的野生型和突变株的表达，比较治疗与未治疗细胞，癌症与正常细胞等。例如，比较培养病原体之间的表达变化鲍曼不动杆菌在乙醇存在下生长的病原体（已知其会增加毒力）揭示了49个不同表达的基因，这些基因属于一系列功能类别[15]. 在这里，我们概述了用于检测RNA-seq中差异表达（DE）的处理管道，并检查了执行分析的可用方法和开源软件工具。我们还强调了几个需要进一步研究的领域。

大多数RNA-seq实验采集纯化RNA的样本，剪切它，将其转换为cDNA并在高通量平台上测序，如Illumina GA/HiSeq、SOLiD或Roche 454[16]. 这个过程产生了数百万个从cDNA片段一端获取的短（25-300 bp）读取。这个过程中的一个常见变体是从每个cDNA片段的两端生成简短的读取，称为“配对-end”读取。这些平台在化学和处理步骤上有很大差异，但无论具体细节如何，原始数据都由一长串短序列和相关质量分数组成；这些构成了本次审查的切入点。

图中概述了用于DE分析的典型RNA-seq管道1首先，读取被映射到基因组或转录组。其次，根据实验的目的，将每个样本的映射读数组装成基因水平、外显子水平或转录物水平的表达摘要。接下来，将汇总的数据与DE的统计测试进行归一化，得出相关基因的排名列表P（P）-值和折叠更改。最后，从这些列表中获得的生物学见解可以通过执行系统生物学方法获得，类似于在微阵列实验中执行的方法。我们在下文中对RNA-seq数据分析的每个步骤的当前可用方法进行了评论。我们没有提供所有可用工具的完整列表，而是侧重于说明方法的常用开源软件示例（表1). 有关RNA-seq分析软件的完整列表，请参阅[17,18].

图1

检测差异表达的RNA-seq分析管道概述.管道中的步骤在红色方框中；管道的方法组成部分显示在蓝色方框和粗体文本中；每个步骤的软件示例和方法（一个非详尽的列表）以蓝色方框中的常规文本显示。表1列出了所示工具和方法的参考。首先，将读数映射到参考基因组或转录组（使用连接文库来映射跨越外显子边界的读数）；映射读取被组装成表达摘要（计数表，显示读取如何位于编码区、外显子、基因或连接处）；对数据进行归一化处理；进行差异表达（DE）的统计测试，生成一系列相关基因P（P）-值和折叠更改。然后可以使用系统生物学方法从这些列表中获得生物学见解。

映射

为了使用RNA-seq数据比较样本之间的表达，有必要将数百万次的简短读取转化为表达的量化。此过程的第一步是读取映射或对齐。最简单的是，映射的任务是找到一个唯一的位置，其中一个简短的读取与引用相同。然而，在现实中，参考从来都不是正在测序的RNA的实际生物来源的完美代表。除了样本特有的属性，如SNP和indels（插入或删除）外，还需要考虑的是，读取是来自剪接转录组而不是基因组。此外，短读有时可以与多个位置完美对齐，并且可能包含必须考虑的排序错误。因此，真正的任务是找到每个短文与参考最匹配的位置，同时考虑到错误和结构变化。

虽然目前正在研究如何将读取与引用对齐，但所有解决方案都必然涉及算法的计算要求和匹配引用时允许的模糊性之间的折衷。几乎所有的短读比对器都使用第一遍“启发式”匹配的策略，该策略可以快速找到可能位置的简化列表，然后通过复杂的“局部比对”算法对所有候选比对进行彻底评估。如果没有此初始启发式搜索来减少潜在对齐位置的数量，则在当前硬件上执行数百万个短读取的本地对齐在计算上是不可能的。

当前对齐器通过使用任意一个哈希表实现快速启发式匹配[19——22]或Burrows-Wheeler变换（BWT）[23——25]. 散列对齐器的优点是易于扩展，以检测读取和引用之间的复杂差异，而代价是不断增加的计算需求。或者，基于BWT的对准器可以非常有效地映射与参考紧密匹配的读数，但一旦考虑到更复杂的不对中，速度就慢得令人望而却步。这些技术的详细解释超出了本次审查的范围，但可以在[23,26——30].

校准器在处理“多重地图”的方式上也有所不同（对多个位置的地图读取效果相同）。大多数对齐器要么放弃多重贴图[25]，随机分配[29]或根据当地覆盖率的估计进行分配[31,32]，尽管也提出了一种包含对齐分数的统计方法[33]. 配对读取减少了多重映射的问题，因为生成短读取的cDNA片段的两端应映射到转录组附近，从而在大多数情况下解决多重映射的模糊性。

当考虑读取基因组DNA时，只需要映射到相关参考基因组。然而，RNA-seq是对转录组的片段进行测序。这种差异可以通过几种方式来处理。鉴于转录组是由基因组“构建”的，最常用的方法（至少在最初）是使用基因组本身作为参考。这样做的好处是容易，不会偏向任何已知的注释。然而，研究表明，跨外显子边界不会映射到该参考。因此，使用基因组作为参考将使具有较少外显子的转录本具有更大的覆盖范围（在相同的真实表达水平上），因为它们将包含较少的外显子连接。较长的读取更有可能跨越外显子边界，从而导致连接读取的比例增加[2].

为了解释连接读取，通常会构建外显子连接库，其中参考序列是使用注释外显子之间的边界构建的[2,32,34,35]. 为了在不依赖现有注释的情况下绘制外显子边界，可以使用数据集本身检测剪接连接从头开始[36——41]. 另一个选项是从头开始转录组的组装，用作参考，使用基因组组装工具[42,43]. 全部从头开始方法可以识别新的转录物，并且可能是没有基因组参考或注释的生物体的唯一选择。然而，从头开始方法是计算密集型的，可能需要较长的成对读取和高水平的覆盖才能可靠地工作。例如，Trapnell等。[11]使用了超过4.3亿个配对读取从头开始组装小鼠成肌细胞转录组，以量化细胞分化过程中的表达。

转录组映射的一种常用方法是逐步增加映射策略的复杂性，以处理未对齐的读取[44]. 例如，在一项调查69个尼日利亚HapMap样本表达变异的大型研究中，Pickrell等。[35]发现对于46 bp的Illumina读数，87%使用MAQ（一种基于散列表的比对器）定位到两个不匹配的参考基因组[23]. 另外7%可以映射到外显子-外显子连接文库，该文库由所有可能的Ensembl外显子组合构建而成。对其余未映射的读取进行了检查，以确定测序器对聚（A）尾进行了错误测序的证据。如果一个读取以至少四个As或T开始或结束，则会修剪这些基数，并将其余读取映射到引用，从而进一步映射0.005%的读取。这个庞大的数据集能够注释100多个新的外显子，并确定了1000多个基因，其中遗传变异影响整体表达水平或剪接。如果没有处理跨越外显子边界读取的方法，这是不可能实现的。

汇总映射读取

在获得尽可能多的读数的基因组位置后，下一项任务是总结和聚合一些有生物学意义的单位的读数，如外显子、转录物或基因。最简单也是最常见的方法是计算与基因外显子重叠的读取次数（例如[32,34,45]). 然而，有很大一部分读码映射到注释外显子以外的基因组区域，即使在注释明确的生物体中，如小鼠和人类。例如，皮克雷尔等。[35]发现约15%的定位阅读位于尼日利亚HapMap样本注释外显子之外，这些外显子阅读更有可能是细胞类型特异性外显子。类似地，图2年显示了RNA结合蛋白39中注释外显子之外发生的转录的示例(RBM39型)LNCaP前列腺癌细胞中的基因。其他正常组织细胞类型的读数更多地局限于已知的外显子，但也显示出已知外显子外转录的证据。

图2

将映射读取汇总为基因级计数.（a）RNA结合蛋白小区域的定位读数39(RBM39型)显示LNCaP前列腺癌细胞的基因[90]UCSC轨道上的人类肝脏和人类睾丸。三行RNA-seq数据（蓝色和黑色图表）显示为“堆积轨迹”，其中每个位置的y轴测量与该位置重叠的映射读取数。还显示了基因组坐标、基因模型（标记为RBM39型; 蓝色方框表示外显子）和脊椎动物的保守性得分。很明显，许多阅读源于没有已知外显子的区域。（b）基因组区域及其可能产生的读码的示意图。读码是根据其来源的基因组特征进行彩色编码的。不同的摘要策略将导致计数表中包含或排除不同的读取集。例如，仅包括来自已知外显子的读取将排除内含子读取（绿色）对结果的贡献。拼接连接被列为一个单独的类别，以强调其赋值中的潜在歧义性（例如应将连接读取赋值给哪个外显子），以及许多读取可能无法映射的可能性，因为它们比连续读取更难映射。CDS，编码序列。

另一种概括是包括整个基因长度的读取，从而合并“内含子”的读取。这将包括摘要中未注释的外显子，并解释注释不足或可变外显子边界。然而，包括内含子也可能捕获重叠的转录物，这些转录物共享一个基因组位置，但来源于不同的基因。还有许多其他可能的变体可用于摘要，例如仅包含映射到编码序列的读取或从从头开始预测外显子[40]. 连接读取也可以添加到基因摘要计数中，或用于建模剪接亚型的丰度[11]. 这些不同的可能性如图所示2亿有了这些选项，摘要的选择就有可能像映射策略的选择一样，或者更大程度上改变每个基因的计数。尽管如此，关于哪种摘要方法最适合DE检测的研究还很少。

规范化

规范化可以准确比较样本之间和样本内的表达式级别[2,32,34]. 研究表明，归一化是从RNA-seq数据分析DE的关键步骤[45——48]. 库内和库间比较的规范化方法不同。

库内标准化允许量化每个基因相对于样本中其他基因的表达水平。由于较长的转录物具有较高的读取计数（在相同的表达水平上），库内标准化的一种常见方法是将汇总计数除以基因的长度[32,34]. 广泛使用的RPKM（外显子模型每千碱基每百万映射读数的读数）解释了样本内比较中的库大小和基因长度效应。为了验证这种方法，Mortazavi等。[32]引入了几个拟南芥RNA进入小鼠组织样本，跨越一系列基因长度和表达水平。这些非天然RNA被称为“尖峰蛋白”，并证明RPKM能够准确比较基因之间的表达水平。然而，研究表明，由于序列内容的原因，表达的转录本上的阅读覆盖可能是不均匀的[49]和RNA制备方法，如随机六聚体启动[50]. 将这种理解融入到内部库规范化方法中可以提高比较表达式级别的能力。使用RNA-seq数据估计样本中转录物的绝对数量是可能的，但它需要RNA标准和其他信息，例如提取RNA的细胞总数和RNA制备产量[32].

当在样本之间测试单个基因的DE时，技术偏差，如基因长度和核苷酸组成，将主要抵消，因为用于汇总的基本序列在样本之间是相同的。然而，样本间归一化对于比较不同库中的计数彼此之间的关系仍然至关重要。最简单和最常用的规范化根据库中的读取总数进行调整[34,51]解释了这样一个事实，即如果一个样本的测序深度更大，每个基因将被分配更多的读取。然而，已经表明，需要更复杂的归一化来解释成分效应[48]或者因为少量高表达基因可以消耗大量的总序列[45]. 为了说明这些特征，可以从数据中估计比例因子，并在测试DE的统计模型中使用[45,46,48]. 缩放因子的优点是保留原始计数数据以供后续分析。或者，分位数归一化和使用匹配幂律分布的方法[52,53]还建议对RNA-seq进行样本间归一化。这两种转换的非线性消除了数据的计数性质，使得如何适当测试DE变得不清楚。到目前为止，分位数标准化似乎并没有像适当的比例因子那样提高DE检测[45]目前尚不清楚幂律分布是否适用于所有数据集[48].

差异表达

DE分析的目的是突出显示在实验条件下大量发生显著变化的基因。通常，这意味着为每个库提取一个汇总计数数据表，并在感兴趣的样本之间执行统计测试。

利用微阵列数据分析差异表达的方法很多。然而，RNA-seq对每个基因进行了离散测量，而微阵列强度具有连续的强度分布。尽管微阵列强度通常作为正态分布的随机变量进行对数变换和分析，但计数数据的变换不能很好地通过连续分布来近似，尤其是在较低的计数范围和小样本中。因此，适用于计数数据的统计模型对于从RNA-seq数据中提取最多信息至关重要。

一般来说，泊松分布是RNA-seq计数数据建模的基础。在一项早期的RNA-seq研究中，使用单一来源的RNA，在Illumina GA测序仪的多个泳道上测序，拟合优度统计表明，大多数基因的计数在泳道上的分布确实是泊松分布[34]. 这已通过技术实验得到独立证实[45]并且可以使用软件工具进行这些分析[54]. 然而，泊松假设并没有很好地捕捉到生物变异性[47,51]. 因此，对具有生物复制的数据集进行基于泊松的分析时，由于低估了抽样误差，容易出现较高的假阳性率[46,47,55]. 尽管RNA-seq平台的背景较低且灵敏度较高，但设计生物复制实验对于确定RNA丰度的变化（普遍适用于采样人群）仍然至关重要。RNA-seq实验的总体设计，包括阻断、随机化和复制的基本考虑，最近已被深入讨论[56].

为了解释生物变异性，基因表达序列分析（SAGE）数据的方法最近被应用于RNA-seq数据[57]. SAGE和RNA-seq数据的主要区别在于数据集的规模。为了解释生物变异性，负二项分布被用作泊松分布的自然延伸，需要估计额外的分散参数。基于阴性肿瘤的计数数据DE分析出现了一些变化，包括常见的离散模型[55]，使用加权似然共享所有基因的信息[51]，均值-方差关系的经验估计[46]和使用等价类的经验贝叶斯实现[58]. 通过广义泊松分布，还提出了泊松模型的扩展，以包括过度分散[59]或两阶段泊松模型，根据数据过度分散的证据，测试两种模式下的差异表达[60]. 同时发现和量化转录物的几种工具[11]或替代亚型表达分析[10]也执行DE分析。然而，值得注意的是，这些方法要么使用泊松分布，要么使用费希尔精确检验，这两种方法都没有明确处理上述生物变化。

许多当前的计数数据DE分析策略仅限于简单的实验设计，如两两或多组比较。据我们所知，在RNA-seq数据的背景下，还没有提出一般方法来分析更复杂的设计，例如配对样本或时间进程实验。在缺乏这种方法的情况下，研究人员已经转换了他们的计数数据，并使用了适合于连续数据的工具[31,47,61]. 广义线性模型为上述计数模型提供了逻辑扩展，需要开发聪明的策略来共享所有基因的信息；软件工具现在提供了这些方法（例如edgeR[57]). 此外，上述方法主要旨在总结注释存在的表达式级别。方法，如最大平均差异测试[62]，最近被提议以非目标方式检测DE。

系统生物学：超越基因列表

在许多情况下，创建DE基因列表并不是分析的最后一步；通过观察多组基因的表达变化，可以进一步了解实验系统的生物学特性。许多专注于基因集测试、网络推理和知识数据库的工具被设计用于分析来自微阵列数据集的DE基因列表[63——65]. 然而，RNA-seq受到微阵列数据中不存在的偏差的影响。例如，基因长度偏差是RNA-seq数据中的一个问题，其中较长的基因具有较高的计数（在相同的表达水平上）[66]. 这使得检测长时间和高表达基因的DE具有更大的统计能力。这些偏差会显著影响下游分析的结果，例如测试DE基因中富集的基因本体（GO）术语[66,67]. 为了进行基因集分析，Bullard等。[45]建议修改DEt吨-通过除以基因长度的平方根来最大限度地减少长度偏差对DE的影响。或者，GO-seq是一种专门针对RNA-seq数据开发的方法，可以将长度或总计数偏差纳入基因集测试[68]. 随着对RNA-seq数据中偏差的理解不断加深，包含这种理解的系统生物学工具对于提取生物学见解至关重要。

将RNA-seq数据的结果与其他生物数据源相结合，以建立一个更完整的基因调控图景有很大的空间[69]. 例如，RNA-seq已与基因分型数据结合使用，以确定导致个体间基因表达变异的遗传位点（表达数量性状位点或eQTL）[35,70]. 此外，将表达数据与转录因子结合、RNA干扰、组蛋白修饰和DNA甲基化信息相结合，有助于更好地理解各种调控机制。最近出现了一些关于这些“综合”分析的报告[71——73]. 例如，李斯特和合著者[71]强调了RNA-seq表达与基因体中CG和非CG甲基化水平的相关性的显著差异。类似地，基于序列的数据集的组合开始提供对表达、组蛋白修饰和DNA甲基化之间单等位基因关联的见解[74].

见解

在这篇综述中，我们概述了将RNA-seq产生的数以百万计的短阅读处理成样本间DE分析的主要步骤。简而言之，该过程是绘制和总结短阅读序列，然后在样本之间进行归一化，并对DE进行统计测试。通过寻找基因组内的表达变化模式，并将RNA-seq数据与其他来源的数据集成，可以获得进一步的生物学见解。

尽管这条管道的许多部分一直是广泛研究的重点，但仍有一些领域可以进一步完善。到目前为止，很少有人研究哪种摘要度量最适合于发现样本之间的DE。此外，还可以扩展现有的DE检测统计方法，以分析更复杂的实验设计。此外，就分析各种研究设计的灵活性、在小型和大型研究中的表现、对测序深度的依赖性以及所施加假设（如均值-方差关系）的准确性而言，目前可用的许多方法的相对优点值得进一步研究。此外，尽管有许多使用RNA-seq检测选择性剪接的例子，但仍有扩大目前检测基因同种型偏好差异的方法的余地[10,11]当生物变异性显著时，可能使用上述基于计数的统计方法。

鉴于产生短读的协议之间存在重大差异，因此正式比较RNA-seq平台和许多数据分析方法的相对优点非常重要。此类调查可能会揭示特定于平台的DE分析方法的优点，也将促进更大的数据集成。由于该领域还相对年轻，我们预计在不久的将来会出现许多用于分析RNA-seq数据的新方法和工具。

方框1：基因表达分析的微阵列和测序比较

现在已经对RNA-seq和微阵列数据进行了一些比较。这些包括测序平台的基本证明[2,31,32]，专门的比较研究[34,75——77]和分析方法开发[10]. 结果一致：测序具有较高的灵敏度和动态范围，加上较低的技术变异。此外，比较突出了微阵列和测序在绝对表达和差异表达测量方面的强烈一致性。尽管如此，微阵列在许多生物环境中询问转录组方面一直并将继续非常成功。例子包括定义乳腺癌亚型的起源细胞[78]并研究进化对基因表达的影响果蝇属[79].

微阵列和测序都有各自的特定偏差，这些偏差可能会影响平台测量DE的能力。众所周知，微阵列探针的交叉杂交以不均匀的方式影响表达测量[80,81]和序列含量影响测得的探针强度[82]. 同时，一些研究观察到RNA-seq数据中存在GC偏差[45]而RNA-seq可能会受到同源序列映射模糊的影响。此外，在较高计数下检测变化的统计能力更高（例如，在无差异的零假设下，200次读取到100次读取的两倍差异比20次读取到10次读取的二倍差异更具统计意义）；这种偏见在RNA-seq中典型地表现为DE和基因长度之间的关联，这种效应在微阵列数据中并不存在[66,68]. 其他研究表明，特定的测序协议会在生成的读取中产生偏差，这可能与序列组成和转录本上的距离有关[49,50,83,84]. 例如，已经发现小RNA的文库准备对观察到的序列集有很大影响[85]. 此外，转录组组装方法必然受到表达水平的影响，因为低水平表达的基因信息较少[11,14]. 许多这些偏见仍在探索中，利用这些知识的聪明统计方法可能能够改进现有方法。

除了RNA-seq的较大动态范围和敏感性外，还有几个其他因素有助于差异表达分析测序的快速应用。首先，微阵列对于许多非模式生物来说根本不可用（例如，Affymetrix为大约30个物种提供微阵列[86]). 相比之下，数千种物种的基因组和序列信息很容易获得[87]. 此外，即使基因组不可用，RNA-seq仍然可以执行，转录组仍然可以被询问（例如，最近的一项研究使用RNA-seq来调查塔斯马尼亚魔鬼面部肿瘤的细胞起源[88]). 其次，测序提供了阵列无法提供的转录特征的前所未有的细节，例如新转录区域、等位基因特异表达、RNA编辑和捕获选择性剪接的综合能力。例如，最近的一项RNA-seq研究[11]能够显示细胞分化过程中的几个亚型转换示例，并且RNA-seq用于显示小鼠大脑中的亲本原表达[5].

当然，排序并非没有挑战。该平台的成本可能会限制一些研究。然而，随着总测序能力和多路复用能力的扩大，每个样本产生足够序列深度的成本很快将与微阵列相当。然而，信息学存储、处理和分析数据的成本是巨大的[89]. 由于数据分析程序相对成熟，计算人员和资源有限的研究人员可能会选择使用微阵列。最后，很明显，测序数据的数据分析方法将在一段时间内继续发展。

作者信息

所有作者对这篇综述的贡献都是一样的。