RNA-seq数据差异表达分析的尺度归一化方法

摘要

基于序列的转录组调查提供的精细细节表明，RNA-seq很可能成为询问稳态RNA的平台。为了发现表达中的重要生物学变化，我们表明正常化仍然是分析中的一个重要步骤。我们概述了一种执行规范化的简单而有效的方法，并显示了在模拟和公共可用数据集中推断差异表达式的显著改进的结果。

背景

转录结构是细胞功能的一个复杂而动态的方面。下一代稳态RNA测序（RNA-seq）提供了细胞内RNA景观的前所未有的细节。不仅可以在没有特定先验知识的情况下查询基因的表达水平，还可以比较样本中基因之间的表达水平。还证明了剪接变体[1,2]和单核苷酸多态性[三]可以通过对转录组进行测序来检测，这为询问等位基因特异性表达和RNA编辑提供了机会。

处理短读排序产生的大量数据的一个重要方面是用于提取和解释信息的处理方法。微阵列数据的经验一再表明，归一化是处理管道的关键组成部分，可以准确估计和检测差异表达（DE）[4]. 标准化的目的是消除数据中出现的系统技术影响，以确保技术偏差对结果的影响最小。然而，生成RNA-seq数据的程序与微阵列数据的程序有根本不同，因此使用的归一化方法并不直接适用。有人建议，“RNA-seq的一个特别强大的优点是，它可以捕获不同组织或条件下的转录组动力学，而无需复杂的数据集标准化”[5]. 我们在此证明，RNA-seq数据分析的现实并非如此简单；正常化通常仍然是一个重要的考虑因素。

目前的RNA-seq分析方法通常通过将给定泳道或文库中的读取次数缩放到实验中所有测序文库的共同值来标准化样本之间的数据。例如，几位作者用包含总读取次数因子的平均值对基因的观察计数进行了建模[6–8]. 这些方法在用于推断差异的分布假设上可能有所不同，但共识是使用模型中的读取总数。同样，对于LONG-SAGE-seq数据，'t Hoen等. [9]使用标度计数的平方根或Vencio的β-二项式模型等. [10]，两者都使用观察到的标记总数。为了正常化，莫拉扎维等. [11]将它们的计数调整为每千基地每百万映射阅读数（RPKM），这表明它“有助于样本内部和样本之间转录水平的透明比较”相比之下，科洛南等. [12]对数转换基因长度标准化计数数据并应用标准微阵列分析技术（分位数标准化和适度t统计量）。苏丹等. [2]通过基因的“虚拟长度”、外显子序列中独特的27米数以及总读取数来规范读取计数。最近，Balwierz等. [13]说明deepCAGE（基因表达的深度测序帽分析）数据遵循近似幂律分布，并提出了一种标准化策略，该策略等同于样本间的读取计数分布。

作为标准化的一种形式，按文库大小进行缩放具有直观的意义，因为预计将样本测序到一半深度将平均获得每个基因对应的读取数的一半。我们认为这适用于RNA群体重复样本之间的正常化。然而，对于许多生物应用来说，库大小缩放过于简单。预期映射到基因的标签数量不仅取决于基因的表达水平和长度，还取决于正在采样的RNA群体的组成。因此，如果一个实验条件中有大量基因是唯一的或高度表达的，那么该样本中剩余基因的测序“不动产”就会减少。如果不进行调整，这种采样伪影可能会迫使DE分析偏向一种实验条件。当前分析方法[6,11]没有明确说明数据的这种比例特性，可能会导致较高的假阳性率和较低的检测真实差异的能力。

这里的基本问题是跨样本比较的适当表达式度量。标准程序是计算每个基因的读数相对于读数总数的比例，并通过转换原始数据或在统计模型中引入常数来比较所有样本的读数。然而，由于不同的实验条件（例如组织）表达不同的RNA库，我们不能总是期望比例可以直接比较。此外，我们认为，在发现有生物学意义的表达变化时，应该认为引导DE调用的欠采样或过采样效应（下文将进一步讨论）是不可取的。下面介绍的归一化方法使用原始数据估计可用于下游统计分析程序的适当比例因子，从而解释RNA-seq数据的采样特性。

结果和讨论

假设情景

估计的标准化因子应确保在两个样本中具有相同表达水平的基因不会被检测为DE。为了进一步强调RNA-seq数据中需要更复杂的标准化程序，请考虑一个简单的思维实验。假设我们有一个比较两个RNA群体a和B的测序实验。在这个假设场景中，假设在B中表达的每个基因在a中表达的转录物数量相同。然而，假设样品A还包含一组数量和表达相同的基因，而这些基因在B中没有表达。因此，样品A的总表达基因是样品B的两倍，也就是说，其RNA产量是样品B大小的两倍。如果没有任何额外的调整，两个样本中表达的基因的平均读取次数将是样本a的一半，因为读取的基因数是样本a基因数的两倍。因此，正确的归一化将样本A调整为2倍。

上述假设示例强调了这样一个概念，即库中特定基因的读取比例取决于整个样本的表达特性，而不仅仅是该基因的表达水平。显然，上面的例子是人为的。然而，在生物甚至技术情况下，需要进行这种规范化。例如，如果RNA样本被污染，表示污染的读数将从真实样本中删除读数，从而减少感兴趣的读数数量，并抵消每个基因的比例。然而，正如我们所证明的，样本之间RNA组成的真正生物差异将是归一化的主要原因。

采样框架

规范化需求的更正式解释使用以下框架Y（Y） _gk公司 作为基因的观察计数克在库中k个从原始阅读中总结，μ _gk公司 作为真实和未知的表达水平（转录本数量），L（左） _克 作为基因的长度克和N个 _k个 作为库的读取总数k个。我们可以对Y（Y） _gk公司 作为：

S公司 _k个 表示样本的总RNA输出。RNA-seq数据分析背后的问题是N个 _k个 已知的是，S公司 _k个 未知，根据RNA组成，不同样品之间可能会有很大差异。如上所述，如果一个群体的总RNA输出量较大，那么相对于另一个样本，RNA-seq实验将对许多基因进行欠采样。

在这个阶段，我们将上述模型中的方差保留为Y（Y） _gk公司 未指定。根据实验情况，泊松似乎适用于技术复制[6,7]负二项式可能适用于从生物复制中观察到的额外变异[14]. 同样值得注意的是，在实践中L（左） _克 通常被吸收到μ _gk公司 参数，并且不在推理过程中使用。然而，已经很好地证明，基因长度偏差在基因表达分析中非常突出[15].

M值归一化方法的修剪平均值

总RNA产量，S公司 _k个 不能直接估计，因为我们不知道每个基因的表达水平和真实长度。然而，两个样本的相对RNA产量，（f） _k个 =S _k个 /秒 _k’ 基本上是一个全球褶皱变化，可以更容易地确定。我们提出了一种经验策略，在假设大多数样本不是DE的情况下，将样本之间的基因总表达水平等同。一种简单而稳健的估算RNA生成比率的方法是使用对数表达比率的加权修剪平均数（修剪平均数M值（TMM））。对于测序数据，我们将基因的对数变换定义为：

和绝对表达水平：

为了稳妥地总结观察到的M值，我们在取加权平均值之前对M值和A值进行了修剪。精度（方差的倒数）权重用于解释这样一个事实，即具有较大读取计数的基因的对数倍变化（实际上是对数相对风险）在对数尺度上的方差较低。有关详细信息，请参见材料和方法。

对于两个样本的比较，只有一个相对比例因子(（f） _k个 )是必需的。它可用于调整两个库大小（将参考除以并将非引用乘以)在统计分析中（例如，Fisher精确检验；有关更多详细信息，请参阅材料和方法）。

通过选择一个样本作为参考并计算每个非参考样本的TMM因子，可以计算多个样本的归一化因子。与双样本比较类似，TMM归一化因子可以构建到用于测试DE的统计模型中。例如，泊松模型会将观察到的库大小修改为有效的库大小，从而调整建模的平均值（例如，在广义线性模型中使用额外的偏移量；有关详细信息，请参见材料和方法）。

肝脏与肾脏数据集

我们将我们的方法应用于公开可用的转录谱数据集，比较肝脏和肾脏RNA来源的几个技术复制[6]. 图1a个显示了在计算读取总数的标准归一化程序后，肾脏样品的两个技术复制品之间的M值分布。这些技术复制品的M值分布集中在零附近。然而，图1亿显示肝脏和肾脏样本之间的对数比显著地向肾脏中更高的表达偏移，即使在考虑读取总数之后也是如此。另一个突出显示的（绿线）是一组家政基因的观察到的M值的分布，显示出显著偏离零。如果按比例缩放到适当标准化RNA-seq数据的读取总数，则预期不会出现这种对数变换偏移。这种偏见的解释很简单。图中的M与A图1c个说明肝脏中存在一组显著的高表达基因（黑箭头）。因此，M值的分布（肝脏到肾脏）呈负向倾斜。由于大量测序都是针对这些肝脏特异性基因的，因此剩余基因的测序较少，因此M值（因此DE值）成比例地扭曲为肾特异性。

图1

RNA-seq数据需要归一化。数据来自[6]比较的对数比（a）技术复制和（b）调整每个样本中读取的总次数后，肝脏和肾脏的表达水平。绿线表示管家基因对数变化的平滑分布。（c）比较肝脏和肾脏的M与A图显示出明显的零偏移。绿点表示545个内务基因，而绿线表示内务基因的中位数对数比率。红线表示估计的TMM归一化因子。橙色斑点涂片突出显示了仅在一个肝或肾组织中观察到的基因。黑色箭头突出显示了一组显著的基因，这些基因在很大程度上归因于对数变换中的总体偏差。

对这对样品应用TMM归一化后，归一化因子为0.68（对数2标度为-0.56；如图中的红线所示1b、c)，反映了大多数肝脏基因的采样不足。TMM因子对于覆盖率较低的数据是稳健的，在这些数据中，可能会有更多的零计数基因（附加文件中的图S1a1)并且对于配平参数的合理值是稳定的（附加文件中的图S1b1). 在DE的统计测试中使用TMM归一化（见材料和方法）会导致肝脏（47%）和肾脏（53%）中相似数量的基因显著增加。相比之下，标准规范化（到最初用于[6])导致大多数DE基因在肾脏中显著较高（77%）。值得注意的是，在TMM归一化后，使用标准归一化方法鉴定为DE的基因中仍检测到不到70%（表1). 此外，我们还发现了大量家政基因的对数变换（来自[16])平均而言，从零开始的偏移非常接近估计的TMM因子，从而使我们的稳健估计程序具有可信度。此外，使用未调整的测试程序，肝脏和肾脏中分别有8%和70%的管家基因显著上调。TMM调整后，DE管家基因的比例分别变为26%和41%，这是一个较低的总数，并且两个组织之间更加对称。当然，假设微阵列数据已经适当归一化（附加文件中的图S21). 综上所述，这些结果表明RNA-seq数据规范化的关键作用。

表1使用不同的归一化方法，以<0.001的错误发现率在肝脏和肾脏之间差异表达的基因数量

其他数据集

在其他RNA-seq数据集中，由RNA组成差异引起的对数变换的全球变化在不同程度上发生。例如，克隆人的M与A情节等. [12]数据集（附加文件中的图S31)给出了在SOLiD™系统上测序的两个样本（类胚体与胚胎干细胞）之间的估计TMM比例因子1.04。该数据集的M与A图也突出了一组有趣的基因，这些基因总体表达较低，但在类胚体中表达较高。这解释了剩余基因对数变换的正向变化。TMM比例因子似乎接近一组大约500个小鼠家政基因（来自[17]). 另一个例子是，李等. [18]使用llumina 1G基因组分析仪的数据集显示出对数变化的总体分布发生了变化，并给出了0.904的TMM比例因子（附加文件中的图S41). 然而，有些基于序列的数据集具有非常相似的RNA输出，可能不需要进行重大调整。例如，来自Kuchenbauer的小型RNA-seq数据等. [19]在log-fold-changes中只显示了适度的偏差（附加文件中的图S51).

尖峰控制有可能用于标准化。在这种情况下，将少量已知量的外来生物体RNA以特定浓度添加到每个样本中。为了使用加标对照进行归一化，在整个实验过程中，加标与样品的浓度比必须保持恒定。在实践中，这很难实现，并且小的变化将导致归一化因子的有偏估计。例如，使用来自Mortazavi的插入DNA等.数据集[11]将导致不切实际的标准化因子估计（附加文件中的图S61). 与微阵列一样使用实验数据仔细估计归一化因子通常更稳健（例如[20]).

模拟研究

为了研究TMM归一化方法的适用范围，我们开发了一个模拟框架来研究RNA组成对RNA-seq数据DE分析的影响。首先，我们模拟来自两个库的数据。我们包括每个样本唯一表达的基因数量参数，以及样本间差异表达基因的比例、大小和方向参数（见材料和方法）。图2a个显示了包含独特基因和DE基因的典型模拟的M与A图。通过模拟不同的总RNA输出，大多数非DE基因具有从零开始偏移的对数变换。在这种情况下，使用TMM归一化来解释潜在的RNA组成，可以减少使用Fisher精确测试的错误检测次数（图2亿). 在广泛的仿真参数范围内多次重复仿真，当比较仿真的真实归一化因子与使用TMM归一化估计的那些因子时，我们发现非常一致（附加文件中的图S71).

图2

仿真表明，TMM归一化是稳健的，优于库大小归一化.（a）模拟结果的示例表明，由于基因在一个样本中唯一表达（橙色点）和不对称DE（蓝色点），因此需要进行标准化。（b）与标准归一化相比，TMM归一化的假阳性率更低。

为了进一步比较TMM归一化与之前在RNA-seq数据DE分析中使用的方法的性能，我们扩展了上述模拟以包括重复测序运行。具体来说，我们比较了三种已发布的方法：经过对数转换的长度规范化计数数据和由Cloonan实现的分位数规范化等. [12]，泊松回归[6]具有库大小、TMM归一化和泊松精确测试[8]具有库大小和TMM规范化。我们没有直接与Balwierz提出的标准化进行比较等. [13]因为肝脏和肾脏数据集似乎不遵循幂律分布，并且具有非常明显的计数分布（附加文件中的图S81). 此外，鉴于我们观察到的RNA组成偏差，尚不清楚等值样本的计数分布是否是最合理的程序。此外，我们没有直接将归一化与虚拟长度进行比较[2]或RPKM[11]归一化，因为没有提到对转换数据的统计分析。然而，我们用M与A图说明，它们的归一化并不能完全消除RNA组成偏差（附加文件中的图S9和S101).

在模拟中，我们使用了基因长度和计数的经验联合分布，因为克隆人等.过程同时需要这两个参数。我们将模拟数据泊松分布以模拟技术复制（附加文件中的图S111). 图3a年显示了两种情况下常见的基因之间的错误发现图，其中我们在第一种情况下引入了10%的唯一组表达，在两倍水平上引入了5%的DE，其中80%在第一种条件下更高。使用为微阵列数据开发的方法学的方法学表现一致较差，正如人们可能预期的那样，因为这些方法的分布假设完全不同。在其余方法（泊松似然比统计、泊松精确统计）中，性能非常相似；同样，TMM规范化对这两者都有显著的改进。

图3

比较几种已发布方法的错误发现图。红线表示长度规范化缓和t统计分析。实线和虚线分别显示了库大小归一化和TMM归一化泊松模型分析。蓝色和黑色线条分别代表LR测试和精确测试。可以看出，使用TMM归一化可以大大降低错误发现率。

结论

TMM归一化是从RNA-seq数据估计相对RNA生成水平的一种简单而有效的方法。TMM方法估计样本之间的标度因子，这些标度因子可纳入当前使用的DE分析统计方法。我们已经证明，在样本之间表达的转录本的潜在分布显著不同的情况下，需要进行标准化。TMM方法背后的假设与微阵列标准化过程中常见的假设类似，如低标准化[21]和分位数归一化[22]. 因此，充分规范化的阵列数据不能显示样本之间不同总RNA输出的影响。本质上，微阵列和TMM归一化都假设两个样本共有的大多数基因没有差异表达。我们的模拟研究表明，TMM方法对一个方向上高达约30%的DE的偏离是稳健的。对于许多应用程序，这一假设不会被违反。

与TMM规范化RNA-seq的一个显著区别是，数据本身不需要修改，不同于微阵列规范化和一些已实施的RNA-seq策略[11,12]. 这里，估计的归一化因子直接用于测试DE的统计模型，同时保留数据的采样属性。由于数据本身没有修改，因此可以用于进一步的应用，例如比较基因之间的表达。

归一化在高通量测序的许多其他应用中至关重要，因为所比较的DNA或RNA群体的组成不同。例如，染色质免疫沉淀（ChIP）之后的下一代测序（ChIP-seq）可能需要进行类似的调整，以比较含有不同结合靶点的样本。有趣的是，PeakSeq方法[23]使用全基因组装箱计数的线性回归来估计两个ChIP种群之间的比例因子，以解释不同的覆盖率。这在原则上与这里提出的类似，但可能不太稳健。我们证明，在许多生物情况下，需要进行成分调整。此外，库大小调整未完全捕获的技术工件可以通过经验调整进行解释。此外，尚不清楚已知浓度的DNA峰值是否能够对归一化因子进行稳健估计。

与以前的高通量技术（如微阵列）类似，归一化是推断样本之间真实表达差异的关键步骤。基因的读取次数不仅取决于基因的表达水平和长度，还取决于其来源的RNA数量。我们提出了一种简单有效的经验方法来规范RNA-seq数据。

材料和方法

TMM规范化详细信息

修剪平均值是删除上下x%数据后的平均值。TMM程序通过对数变换进行双重调整（样品k个相对于样品第页对于基因克)和绝对强度(A类 _克 ). 默认情况下，我们会修剪M（M） _克 值增加30%A类 _克 值减少5%，但这些设置可以根据给定的实验进行定制。该软件还允许用户设置a值的下限，例如Cloonan等.dataset（附加文件中的图S11). 修剪后，我们取加权平均值M（M） _克 ，权重是近似渐近方差的倒数（使用delta方法计算[24]). 具体来说，样本的归一化因子k个使用参考样品第页计算如下：

以下情况Y（Y） _gk公司 =0或Y（Y） _克 ＝0在该计算之前被修剪，因为不能计算对数折叠变化；G公司*表示一组有效的基因M（M） _克 和A类 _克 值，而不是使用上述百分比进行修剪。应该清楚的是.

如图所示2a个表明，总计数越高，M值的方差越低。在一个库中，计数向量是多项式分布的，任何单个基因都是二项分布的，具有给定的库大小和比例。使用增量方法，可以计算M（M） _克 ，这与对数相对风险通常所做的一样，其倒数用于加权平均值。

我们在一系列仿真参数上比较了加权和未加权的修剪平均值以及另一种稳健估计（稳健线性模型），如附加文件中的图S4所示1.

管家基因

人类管家基因，如[16]，已从下载[25]并使用生物导体与Ensemble基因标识符匹配[26]biomeRt包[27]. 类似地，根据deJonge的计算，小鼠家政基因被认为是变异系数最低的大约500个基因等. [17].

统计测试

对于两个库的比较，我们使用CRAN statmod包中的sage.test函数[28]计算Fisher精确P（P）-每个基因的值。为了应用TMM规范化，我们将原始库大小替换为“有效”库大小。对于两个库，通过将估计的归一化因子的平方根与原始库大小相乘/相除来计算有效库大小。

为了与技术复制品进行比较，我们遵循了Marioni中使用的分析程序等.研究[6]. 简单地说，根据以下公式，假设映射到基因的计数是泊松分布的：

哪里表示基因读取总数的分数克在实验条件下z（z） _k个 他们的分析使用偏移量来解释库大小，并使用似然比（LR）统计来测试库之间的表达差异（即H₀:μ_克1=μ_克2). 为了使用TMM归一化，我们使用估计的归一化因子来增加原始偏移量。然后使用相同的LR测试框架计算P（P）-组织之间的DE值。我们修改了这一分析，以使用精确的泊松检验来测试两个复制组之间的差异。该策略在原则上与Fisher的精确测试类似：以总计数为条件，我们计算出观察到组计数的概率为或比我们实际观察到的更极端。总数和组总数均为泊松分布。

我们重新实现了Cloonan的方法等. [12]使用自定义R分析模拟数据[29]脚本。

模拟详细信息

设置模拟是为了从给定的读取计数经验分布（即，从观察到的Y（Y） _克 ). 平均值由采样读取计数除以总和计算得出S公司 _k个 并乘以指定的库大小N个 _k个 （根据模型）。然后从给定平均值的泊松分布中随机抽样模拟数据。我们有参数来指定两个库中通用的基因数量以及每个样本中唯一的基因数量。附加参数指定DE的数量、方向和大小以及测序深度（即读取总数的范围）。由于我们已经插入了已知的差异表达基因，我们可以根据各种统计数据对基因进行排序，并将错误发现的数量作为排序的函数绘制出来。附加文件中的表S11给出了用于图中所示仿真的参数设置2和三.

软件

实现我们方法的软件在edgeR包中发布[30]在Bioconductor版本2.5中[26]和可从[31]. 我们的分析脚本和数据，包括模拟框架，可从[32].