From RNA-seq reads to differenti

image.png

Alicia Oshlack, Mark D Robinson, Matthew D Young.
From RNA-seq reads to differential expression results.Genome Biology 2010, 11:220 http://genomebiology.com/2010/11/12/220

这篇评论文章摘要只有一句话:有很多可用的方法和工具进行预处理高通量RNA-seq数据和检测差异表达。

高通量测序技术现已普遍用于生物学中。这些技术产生了数以百万计的短序列reads，常规地应用于基因组、表观基因组和转录组。对样品中稳定状态的RNA 进行测序，称为RNA-seq，从先前技术的很多限制中解放出来，如对先验的物种知识的依赖，这是微阵列和PCR所须的。此外，RNA-seq允许我们阐 明先前难以见到的转录组复杂性，如等位点特异的表达和新启动子、亚型。但是，产生的数据集是大而复杂的，解释不是直接的。正如任何高通量技术一样，分析方 法对解释数据是至关重要的，而RNA-seq分析过程一直在不断演变。因此，是时候回顾当前可用的分析方法和评论问来研究方向了。

理解RNA-seq数据依赖于感兴趣的科学问题。例如，决定等位点表达的差异需要精确确定转录的SNPs的广泛存在。另一方面融合基因或癌症样品中的畸变可以通过寻找RNA-seq数据中的新转录本来检测。过去一年（即2009年），一些方法涌现出来，用RNA-seq数据进行丰度估计，可变剪接、RNA编辑和新转录本的检测。然而，很多生物学研究的基本对象是样品间的基因表达谱。因此，在本评论中我们聚焦于可以检测样品间基因表达水平差异的可用方法。这种分析与控制实验尤其相关，如比较同一组织的野生型和突变株的表达，比较处理与未处理的细胞，癌症的和正常的细胞等等。我们在此列出用于检测RNA-seq差异表达的处理流程，并检查可以执行该分析的可用方法和开源软件。我们还突出了需要进一步研究的一些方面。

多数RNA-seq实验取一个纯化了的RNA样品，切碎，转换成cDNA并在高通量平台如Illumina GA/HiSeq、SOLiD或Roche 454上测序。该过程产生了来自cDNA片段的一端的、数以百万计的reads（25~300bp）。该过程一个常用的变式是生成双末端reads，即paired-end reads。各平台在化学和处理步骤上有本质不同，但忽略精确的细节后，原始数据都是由带有质量值的短序列的一个长列表组成；这就形成了本评论的切入点。

图1列出了差异表达分析的典型RNA-seq流程概览。首先，reads映射到基因组或转录组。其次，每个样品映射的reads依实验目的而组装成基因水 平、外显子水平或转录本水平的表达概括。接下来，汇总的数据进行归一化以与差异表达的统计检验相协调，产生了一个带有P-value和倍数变化的、排好序 的基因列表。最后，执行系统生物学方法从这些列表中获得生物学的见解，就像在微阵列上进行的那样。我们批判了下列目前可用的RNA-seq数据分析方法的 每一步。我们聚焦于普遍可用的开源软件而不是提供一个所有工具的完整列表。

映射

为了用RNA-seq数据比较样品间表达水平，必须把短reads转换成表达定量。这个过程的第一步就是read映射或比对。最简单的，映射工作就是找到 短read与参考序列已知的唯一位置。然而，真实情形中参考序列从来不是所测序RNA的实际生物源的完美表示。除了样品特异的属性如SNPs和 indels之外，还要考虑来自剪接过的转录组而非基因组的reads。而且，短read有时完美地比对到多个位置，也可能包含不得不考虑的测序错误。因 此，真正的任务是找到短read最佳匹配到参考序列的位置，其中允许错误和结构变异。

尽管对如何最佳比对reads到参考序列的研究还在进行，但是所有的解决办法都涉及在算法的计算需求和允许匹配参考序列的模糊性之间一定的妥协。几乎所有 的短read比对器都采用了首先通过启发式匹配的策略，这迅速找打了可能位置的一个简化列表，接着对候选位置进行全面评估，通过一个复杂的局部比对算法。 如果不做预先的启发式搜索来约减潜在的比对位置数，那么在目前的硬件上执行百万级短reads的局部比对会是计算上不可能的。

目前的比对器能用hash表或Burrows Wheeler变换（BWT）进行快速启发式匹配。hash表比对器对于检测read和参考序列的复杂差异有易于扩展的优点，在不断增加计算需求的代价 下。而BWT比对器可以很有效地映射很接近匹配参考序列的reads，但是一旦考虑更复杂的错配就会大幅度慢下来。这些技术的详细说明可参考文献 23,26-30.

比对器在怎样处理多映射方面也很不一样。多数比对器要么忽略多映射、随机定位它们，要么基于局部覆盖度的估计来定位，尽管结合比对分数的方法也已经提出。PE reads减少了多映射问题，使得多映射的模糊性在大多数情况下都可以解决。

当考虑reads来自基因组DNA时，所有要做的就是映射到一个相关的参考基因组上。但是，RNA-seq是测序转录组片段。这个差异可以用几种方式处 理。既然转录组是建立在基因组之上，那么最常用的（至少是最初的）方法就是用基因组自身做参考序列。这有容易而不偏向任何已知注释的好处。但是跨外显子边 界的reads不会映射到参考序列上。因此，用基因组做参考序列会给出较少外显子的转录本以更高的覆盖度。越长的reads越可能跨外显子边界，因此引起 接合reads比例增加。

为考虑接合reads，通常的实践是建立外显子接合位点库，其中用注视了的外显子的边界构建了参考序列。为了不依赖现存注释的跨外显子边界，可用数据集自 身来从头检测剪接位点。另一个选择是转录组从头组装。所有的从头方法都能鉴定新转录本，并且对没有参考基因组或注释的物种来说是唯一的选项。但是从头方法 是计算密集的，需要长PE reads和高覆盖度以可靠地进行计算。

常用的转录组映射方法是逐渐增加映射策略的复杂性以处理未比对上的reads。

汇总映射的reads

已经尽可能多的获得了reads的基因组位置，下一步任务就是在一些生物学意义单位上汇总和合计这些reads，如外显子、转录本、基因等水平。最简单的 最常用的方法是计数与基因的外显子重叠的reads。但是有相当部分reads映射到基因组上已注释外显子以外的区域，即使是良好注释的物种，如小鼠、 人。

一个可选的汇总是包含沿基因全长的reads，从而结合了内含子reads。这就在汇总中包含了未注释的外显子并考虑了注释不太好或可变的外显子边界。但 是，包含内含子也会捕获到重叠转录本——它们共享一个基因组位置但是源于不同基因。还有其他很多可能的变体用于汇总，如只包含映射到编码序列的reads 或者只汇总从头预测的外显子的reads。接合reads也可添加到基因汇总计数中或用于对剪接亚型的丰度进行建模。这些不同的可能性在图2b中图示说 明。在这些选项下，汇总的选择可能大幅改变每个基因的reads计数，甚至比映射策略的选择影响要大。尽管如此，很少有研究实现了哪种汇总方法是最适合差 异表达检测的。

归一化

归一化使得可以比较样品间和样品内的表达水平。已经证明，归一化是RNA-seq数据的差异表达分析的一个关键步骤。文库内和文库间比较的归一化方法是不同的。

文库内归一化允许定量每个基因相对于样品内其他基因的表达水平。因为越长的转录本有越高的reads计数，文库内归一化的常用方法是用基因长度去除汇总计 数。广泛使用的RPKM在样品内比较中同时解释了文库大小和基因长度的影响。为了验证该方法，Mortazavi等引入了一些阿拉伯芥RNAs到小鼠的组 织样品中，跨过一系列基因长度和表达水平。这些非天然的RNAs称为spike-ins，说明了RPKM给出了基因间表达水平的精确比较。然而，已经证明 表达的转录本的read覆盖深度是不一致的，因为序列内容和RNA制备方法，如随机六聚体引发。把这些认识结合到文库内归一化方法中可能会改进比较表达水 平的能力。使用RNA-seq数据来估计样品中转录本的绝对数目也是可能的，但是这需要RNA标准品和额外信息，如总细胞数和RNA制备产出率。

在样本间检验单个基因的差异表达时，技术偏倚如基因长度与核酸组成，大部分会抵消，因为用于汇总的基础序列在样本之间是相同的。然而，样本间归一化对于相 对不同的文库的比较计数来说仍是必要的。最简单最常用的归一化通过文库的总reads进行调整，考虑了测序深度的影响。但是已经证明需要更聪明的归一化来 考虑组成的影响，或这说事实上小部分高表达基因会占总序列数的相当部分。为了对这些特征进行说明，可以从数据中估计出尺度因子，用于差异表达检验的统计模 型。对于后续分析来说，尺度因子比原始计数有优势。另一方面，分位数归一化和一种用匹配指数律分布的方法也被提出用于RNA-seq的样本间归一化。这些 变换的非线性去除了数据的计数本质，使得不清楚怎样合适地进行差异表达检验。目前，分位数归一化似乎并未改善差异表达检测到合适的尺度因子那样的程度，也 不清楚指数律分布应用于所有数据集的情况。

差异表达

差异表达分析的目标是突出在不同实验条件下丰度显著变化的基因。一般地，这意味着得到每个文库的汇总计数数据表并进行感兴趣样本间的统计检验。

很多方法开发出来以进行微阵列数据的差异表达分析。然而，RNA-seq给出每个基因的离散度量而微阵列的强度给出了一个连续的强度分布。尽管微阵列的强 度通常是对数变换过的，而作为正态分布的随机变量进行分析，但是计数数据的转换并不能用连续分布很好的逼近，特别是对低计数范围和小样本。因此，合适于计 数数据的统计模型对于抽取RNA-seq数据的大部分信息是很重要的。

一般地，Poisson分布构成了对RNA-seq数据进行建模的基础。在任何使用单RNA源的RNA-seq研究中，在一个Illumina GA测序仪的多个lane上进行测序，拟合优度检验表明多数基因在lane之间的分布事实上是Poisson分布的。这一点被独立的技术试验验证过，而且 已有可用的软件工具执行该分析。但是Poisson假设并未很好滴捕捉到生物学变异。因此对于有生物学重复的数据集进行基于Poisson的分析将倾向于 高的假阳性率，源于低估了取样误差。尽管RNA-seq平台有低背景、高敏感性，但是带有生物学重复的试验设计对于将RNA丰度变化推广到取样群体中仍是 至关重要的。一般，RNA-seq试验设计，包括分组、随机化和重复都已经深入讨论过。

为了解释生物学变异性，为SAGE数据而开发的方法最近被用于RNA-seq数据。两者的主要差异在于数据集的规模。为了解释生物学变异性，使用负二项分 布作为Poisson分布的自然推广，要估计一个额外的散度参数。一些基于负二项分布的计数数据差异表达分析变体涌现出来，包括普通散度模型，用加权似然 率共享所有基因的信息，均值方差关系的经验估计和使用等价类的经验Bayes实现。Poisson模型的包含散度的扩展也已提出，如广义Poisson分 布或两步Poisson模型（在取决于数据中散度证据的两个模式下检验差异表达）。一些实时转录本发现和定量或可变亚型分析的工具也执行差异表达分析。但 是值得注意的是这些方法要么用Poisson分布要么用Fisher 精确检验，显然都不能处理上面讨论的生物学变异。

当前的很多计数数据差异表达分析策略中有很多都限于简单实验设计，如成对或多组比较。据我们所知，还没有针对更复杂设计的分析的通用方法提出来，如配对样 本会时间过程实验，在RNA-seq数据的语境下。缺少这样的方法，研究者就把他们的计数数据转换成有合适工具的连续数据。一般线性模型提供了对上述计数 数据的逻辑推广，而也需要开发更聪明的策略来共享所有基因的信息；目前软件工具提供了这些方法。进一步，上面讨论的方法主要目的在于汇总已有注释的表达水 平。以无目标方式检测差异表达的方法近来也提出来了，如极大均值差异检验。

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

在很多情形中，建立差异表达基因列表并非分析的最终步骤；可以通过寻找基因集的表达变化获得对试验系统的深入生物学见解。很多聚焦于基因集检验、网络推断 和知识库的工具为分析微阵列数据集的差异表达基因而设计出来。然而，RNA-seq受到微阵列数据所没有的一些偏倚所影响。例如，基因长度偏倚是RNA- seq数据的一个问题，其中越长的基因有越高的计数。这导致了对长的高表达基因来说，差异表达检测有更高的检定力。这些偏倚极大地影响下游分析结果，如 GO富集。为了能进行基因集分析，Bullard等建议通过除以基因长度的平方根来修改差异表达t-statictic以最小化差异表达中长度偏倚的影 响。另一方面，GOseq特别为RNA-seq数据开发的工具，可以把长度或总计数偏倚合并到基因集检验中。随着对RNA-seq数据偏倚的认识深化，结 合了这种认识的系统生物学工具对于提取出生物学见解是至关重要的。

对于集成RNA-seq数据结果和其他生物学数据源已建立更完整的基因调控图，有着广泛的理解。例如，RNA-seq与基因分型数据结合以鉴定解释个体间 基因表达变异的遗传基因座（表达数量性状位点，eQTLs）。而且，整合表达数据和转录因子结合、RNA干扰、组蛋白修饰以及DNA甲基化信息具有更好理 解各种调控机制的潜力。这种整合性分析的一些报道近来也出现了。例如，Lister等突出了基因体中CG和非CG甲基化水平与RNA-seq表达水平的显 著差异。类似地，测序数据集的组合正开始提供单等位基因与表达、组蛋白修饰和DNA甲基化之间关联性的见解。

Outlook

本评论中，我们列出了处理RNA-seq短reads以进行样本间差异表达分析的主要步骤。简言之，过程就是，映射并汇总短reads序列，然后样本间归 一化并执行差异表达统计检验。进一步的生物学见解可以通过寻找基因集内表达变化模式和整合RNA-seq数据和其它来源的数据来获得。

尽管这个流程的很多部分都是扩展研究的焦点，但仍有些领域存在进一步细化的可能。目前，很少有工作在研究那种汇总度量是最适合寻找样本间差异表达基因的。 为了进行更复杂试验设计的分析，还有扩展现有差异表达检测统计方法的余地。而且，现有的很多方法的相对优点应在进一步的研究中经受考验，依照其分析各种研 究设计的灵活性，其在大大小小的研究中的性能，对测序深度的依赖和强加的假设（如均值方差关系）的准确性。进一步，尽管有很多用RNA-seq进行可变剪 接检测的例子，但是仍有必要扩展当前方法以在生物学变异占主导时检测基因亚型偏好的差异，也许使用上述基于计数的方法。

由于产生短reads的实验协议之间有本质不同，正式比较RNA-seq平台以及很多数据分析方法的相对优点将是很重要的。这样的研究会揭示平台特异的差 异表达分析的好处并会促进更好的数据整合。该领域仍然相对年轻，我们希望在未来有更多的RNA-seq数据分析新方法和工具涌现出来。