mRNA-Seq实验中归一化和差异表达的统计方法评估

摘要

背景

高通量测序技术，如Illumina基因组分析仪，是研究广泛生物和医学问题的强大新工具。统计和计算方法是从测序器生成的大量复杂数据集得出有意义和准确结论的关键。我们对Illumina转录组测序（mRNA-Seq）数据的归一化和差异表达（DE）分析的统计方法进行了详细评估。

结果

我们比较了两种生物样本中检测显著DE基因的统计方法，发现测试统计数据处理低计数基因的方式存在显著差异。我们评估了测序平台的特征对DE结果的影响，例如，不同的基因长度、碱基校准方法（有和没有φX控制通道）以及流细胞/文库制备效果。我们研究了读取计数归一化方法对DE结果的影响，并表明通过总车道数（例如RPKM）进行缩放的标准方法可能会对DE估计值产生偏差。我们提出了更通用的基于分位数的归一化过程，并证明了DE检测的改进。

结论

我们的结果对mRNA-Seq实验的设计和分析具有重要的实用性和方法学意义。他们强调了归一化和DE推断的适当统计方法的重要性，以说明测序平台可能影响结果准确性的特征。它们还揭示了在发展mRNA-Seq的统计和计算方法方面需要进一步研究。

背景

过去十年来，微阵列一直是高通量基因表达研究的首选检测方法。全基因组测序的效率、质量和成本的最新改进促使生物学家迅速放弃微阵列，转而使用超高通量测序，即第二代或下一代测序：例如，Applied Biosystems的SOLiD、Helicos BioSciences的HeliScope、Illumina的基因组分析仪、，罗氏公司的454个生命科学测序系统。这些高通量测序技术已经应用于监测全基因组转录水平（mRNA-Seq）、DNA-蛋白质相互作用（ChIP-Seq）、染色质结构和DNA甲基化[1——9].

我们使用微阵列质量控制（MAQC）项目的参考样本，评估用mRNA-Seq推断差异表达（DE）的统计方法[10]. 利用大约1000个基因的相应定量实时聚合酶链式反应（qRT-PCR）数据，我们比较了不同的标准化和DE程序，并评估了与测序技术相关的可能偏差。对于在被比较的两个样本中表达良好的基因，所检查的测试（Fisher精确测试和基于GLM的测试）无法区分。然而，即使只有一个样本产生低读取计数（例如≤10），它们给出可靠DE估计的能力也存在很大差异。Illumina平台的一个固有偏见是较长基因的优先排序[11]. 考虑到这里的测试，更长的基因更有可能被宣布为DE。我们证明，通过基因长度加权DE统计可以减轻这种影响。

虽然由于流细胞/库制备的差异可以观察到小的“麻烦”技术效果，但我们表明，这些不会对MAQC数据集的差异表达调用产生实质性影响。我们还发现，不在每个流动池中使用标准φX控制通道，如Illumina推荐的校准程序，不会对DE检测产生负面影响。此外，没有phi X通道的自动校准增加了映射读取的数量和质量。在这方面，使用φX车道没有明显的好处；取消这样的控制车道将导致更平衡和更具成本效益的设计。

我们证明，对DE检测的最大影响是归一化程序的选择。由于不同的车道具有不同的总读取计数，即。，测序深度，通常的方法是通过总车道数来衡量每条车道内的基因数：例如，现在的标准是外显子模型每千基每百万映射读数（RPKM）[7]或超几何模型[6]. 我们表明，这种形式的全球标准化受到相对较小比例的高表达基因的严重影响，因此，如果这些少数基因在比较条件下差异表达，则可以给出DE的有偏估计。我们提出了另一种更稳健的基于分位数的归一化过程，可以在不引入额外噪声的情况下消除偏差。

方法

MAQC数据集

本文考虑了与MicroArray质量控制项目相关的两个mRNA-Seq数据集[10]并使用Illumina的基因组分析仪II高通量测序系统获得[12]. 实验分析了两个生物样本：Ambion的人脑参考RNA和Stratagene的人类通用参考RNA，这里分别称为brain和UHR。

在第一个实验（MAQC-2）中，对两种类型的生物样品（脑和UHR）进行了分析，每种样品使用分布在两个流动细胞上的七条通道。两种生物样品各使用一种文库制剂。因此，生物效应是困惑的由于文库准备效应，即大脑和UHR之间mRNA-Seq测量值的一些差异可能仅是由于实验伪影。在第二个实验（MAQC-3）中，使用来自两个流动细胞的14条通道分析了四种不同的UHR文库制剂；每个文库制剂仅在一个流动细胞上进行分析。因此，库准备效果是嵌套的流细胞内效应和流细胞之间的差异与库准备效应相混淆（参见[附加文件1：实验设计的补充图S1]）。MAQC-2和MAQC-3实验的测序读数已保存到短读档案馆，登记号为SRA010153.1。

作为最初MAQC项目的一部分，还选择了大约1000个基因进行qRT-PCR检测[13]. 我们使用这些qRT-PCR数据作为黄金标准来衡量mRNA-Seq测定的基因表达值。此外，还进行了大量的微阵列实验。我们将mRNA-Seq测量值与来自一组Affymetrix人类基因组U133 Plus 2.0阵列（GSE5350，样本AFX_1_[a-B]）的测量值进行了比较[1——5]；参见[附加文件2：补充第S1.2和S1.2节]，了解qRT-PCR和阵列分析的详细信息。

Illumina测序平台概述

我们对Illumina mRNA-Seq实验中涉及的步骤进行了简要的非技术性概述[12]. 对感兴趣的样本进行文库准备，这是一系列步骤，将输入的RNA转化为DNA的小片段，然后通过Illumina机器进行测序。具体来说，从任何总RNA样本开始，Illumina的mRNA-Seq文库制备方案包括poly-A RNA分离、RNA片段化、使用随机引物反转录到cDNA、适配器连接、凝胶大小选择和PCR富集[[14]，图六]。生成的cDNA图书馆被放在八个中的一个车道的流动细胞单个cDNA片段附着在车道表面，随后经过扩增步骤，从而转化为集群双链DNA。然后将流细胞放入测序机中，在测序机上对每个簇进行并行测序。具体来说，在每个周期添加四个荧光标记的核苷酸，并记录每个簇发出的信号。对于每个流动细胞，在给定的循环数下重复此过程，例如，MAQC实验中的35个循环。然后将荧光强度转换为棒球循环数决定了读取; 集群的数量决定了读取的数量。

测序数据的预处理

对于两个MAQC实验，使用Illumina的标准基因组分析仪1.3版预处理管道获得了35个碱基-对-长读数[12,15]. 我们使用Bowtie将读取结果映射到基因组（GRCh37组装）[16].

Illumina的默认基本呼叫算法Bustard可以通过两种方式进行校准。Illumina推荐的方法是为每个流细胞保留一条通道用于DNA测序（通常为phi X DNA），并使用此控制通道中的数据来确定其他七条通道的基本细胞数和质量分数[[15]，补充信息，第7页]。Bustard也可以使用自动校准方法运行，该方法以类似于phred基本调用的方式对基本调用进行评分[17]并且不需要每个流单元都有一条控制车道。在两个MAQC实验中，每个流细胞的一条通道被保留用于对phi X基因组DNA进行测序。对于一个实验（MAQC-2），我们获得了自动校准和phi X校准的读数。

除了讨论基数校准方法的影响的部分外，我们重点关注phi X校准、纯度过滤读取，该读取与基因组唯一对应，最多有两个不匹配。限制阅读映射到基因组意味着外显子-外显子连接阅读被排除在外（约10%的阅读）。此外，库准备协议不允许考虑特定于股的计数，因此映射到正向股和反向股的读取被合并。

单交叉基因的定义

在我们对DE的评估中，我们关注的是基因的整体表达，而不是异构体特异性表达。没有标准的技术来总结具有几种亚型的基因的表达水平（例如，参见[6]和[7]对于不同的方法）。对于给定的基因，我们首先定义一个组成型外显子作为属于基因每个亚型的一组连续外显子碱基（即外显子的一部分或整个外显子）。然后我们定义一个单交（UI）基因作为一个复合基因级感兴趣区域，由组成外显子的联合组成，这些外显子与其他基因的编码外显子不重叠（基于Ensembl，版本55；参见[附加文件2：补充章节S2]）。我们保留了所有通过染色体1-22、X和Y识别的基因。除了包括蛋白质编码基因外，UI基因还代表了许多其他类别的集合注释，例如假基因和小RNA。

规范化

为了推导基因表达测量值并在（几组）车道之间比较这些测量值，首先需要对读取计数进行标准化，以调整不同的车道测序深度和潜在的其他技术影响。这里考虑的所有归一化方法中只有一种是全球的程序，在某种意义上，只有一个因素d日 _我 用于缩放计数（每行）。

我们评估了三种类型的全局归一化：（1）总车道数，如RPKM[7]，（2）预期在生物条件下持续表达的“内务管理”基因的per-lane计数，例如，POLR2A，（3）至少在一条通道中读取的基因的per-lane上四分位基因计数。为了使规范化表达式度量具有可比性，缩放因子本身进行了缩放，以便它们在所有车道上的总和等于所有14条车道上的总计数之和（参见[附加文件2：补充章节S4]）。

表达式量化问题可以用广义线性模型（GLM）来描述，

（1）

其中，读取计数的预期值的自然对数X（X） _{i、 j个} 对于j个th基因在我thlane被建模为基因表达水平的线性函数λ_{一(我)，j个}对于生物条件一(我)在巷道中进行化验我加上偏移量（对数d日 _我 )以及可能的其他技术效果(θ _{i、 j个} ).

最后，我们提出了一个分位数归一化过程，该过程受以下微阵列归一化方法的启发[18]及其在R包aroma.light中的实现。具体来说，对于每条车道，读取计数的分布与根据分拣车道的中值计数定义的参考分布相匹配。对规范化数据进行四舍五入，以生成可与下文描述的DE统计数据一起使用的整数值。

差异表达

我们比较了三种推断DE的方法，每种方法都为每个基因生成一个测试统计量：Fisher精确测试统计量，基于广义线性模型的似然比统计量，如方程式所示(1)、和t吨-基于相同GLM的估计参数的统计。两个不同的t吨-对统计进行评估，使用不同的技术估计估计参数的方差。我们还通过添加参数来评估流细胞效应的影响θ _{i、 j个} 在GLM中或通过Mantel-Haenszel测试，Fisher精确测试的扩展（参见[附加文件2：补充部分S5]）。所有考虑的DE统计数据都可以通过偏移来适应全局标准化d日 _我 。对于基于GLM的统计，偏移量按等式处理(1). Fisher精确检验和Mantel-Haenszel检验比较了j个th基因与d日.

似然比统计是最通用的，因为它们可以用于比较任何数量的生物样本类型，并调整一般实验效果以及样本协变量，例如RNA质量。这个t吨-统计数据只适用于测试两组之间的差异。这个t吨-统计和似然比统计基于同一GLM的最大似然估计，但在某些情况下具有不同的性能。在渐近理论下导出了所有基于GLM统计量的分布性质；因此，对于少量输入样本或低计数，它们的行为可能较差（尽管这不是我们遇到的情况）。相比之下，Fisher的精确检验没有对样本大小进行假设；然而，它只针对全球实验效应进行调整，甚至曼特尔-海恩斯泽尔扩展也只允许单一基因水平的实验效应。

似然比统计已用于[6]对于只有全局车道效应的特殊情况（即。，θ _{i、 j个} =方程式中的0(1)); 这些作者还提到了应用反正弦根变换来稳定每个泳道内的每个基因读取比例的方差。泊松参数的Gamma先验贝叶斯统计结果与上述基于GLM的测试统计结果类似[19]. 近期mRNA-Seq文献中考虑的其他测试统计数据包括t吨-基于贝塔二项分布的平方根转换标准误差统计和贝叶斯统计[三].

使用qRT-PCR金标准的接收机操作员特性曲线

的qRT-PCR数据[13]被用作金标准，以确定各种信使核糖核酸序列和微阵列DE方法的“真实”差异表达，并导出受体-操作员特征（ROC）曲线。UHR与大脑表达对数变化的qRT-PCR估计值是UHR和大脑在不同复制品中平均表达测量值的差异（参见[附加文件2：补充章节S6]）。

我们将qRT-PCR分析的基因分为“non-DE”、“DE”和“no-call”三组，根据它们的绝对表达log-fold-change是否小于一，大于b条，或在间隔内[一,b条]分别是。在确定真/假阳性/阴性时，我们忽略了“no-call”基因。当测序（或微阵列）平台不仅正确地声明了一个基因DE，而且在DE的方向上与qRT-PCR一致时，就会报告真阳性（TP）。根据qRT-PCR，真阳性率（TPR）定义为TPs总数除以DE基因总数；像往常一样计算假阳性率（FPR）。请参阅表1以获取摘要。

表1真阳性率和假阳性率的定义。调用真/假阳性和阴性的规则概要，其中考虑了差异表达方向的符号：“+”表示UHR中的过度表达，“-”表示大脑中的过度表示。

软件

为了便于mRNA-Seq数据的分析和可视化，我们开发了两个R/Bioconductor软件包，Genometer和GenomeGraphs[20]. 这两个包都可以从生物导体项目中获得，http://bioconductor.org/packages/release/bioc/html/Genominator.html和http://bioconductor.org/packages/release/bioc/html/GenomeGraphs.html分别是。

结果和讨论

mRNA-Seq差异表达统计的比较

差异表达基因的列表通常是通过计算每个基因的测试统计来生成的，该测试统计比较了两种生物样本之间的表达水平，并根据第页-评估观察差异的统计显著性的值。

我们评估了差异表达的各种统计数据（见上文方法中的描述），发现测试统计数据之间的主要差异在于它们处理低计数的能力，这在研究mRNA-Seq背景下的差异表达时是一个非常重要的问题。当两个样本的读数都为零时，显然无法说明微分表达式。当一个基因对一个样本的读数为零，而对另一个样本的读数为合理数字时，就会出现更相关的零计数或低计数情况。大约700个基因（约1.8%）在大脑或UHR中为零，在其他组织中为10个或更多。据推测，这代表了一种有趣的生物现象，根据测序，一个组织中的基因完全不受压。

对于任一样本中计数为零的基因t吨-统计失败：估计的标准误差变得非常大（在delta方法中为无穷大t吨-统计）和标称值第页-值聚集在一个样本周围，而与另一个样本中的读取次数无关。然而，对于Fisher精确检验和基于GLM的似然比检验，我们可以看到第页-值。对于两个样本中计数合理的基因，测试统计量的选择对标称值几乎没有影响第页-值（[附加文件1：补充图S2和S3]）。因为它们不能稳定地处理低计数基因t吨-统计数据未能检测到许多“简单”的DE病例（即两种情况下表达差异较大的基因），因此敏感性很低。The poor performance of Thet吨-统计数据反映在ROC曲线中（图1). 去除两个样本中读数均小于20的基因完全解释了t吨-统计数据和各种DE统计数据的等效ROC结果，所有这些都得到了显著改进（图1).

图1

差异表达统计比较：ROC曲线（a）所有DE统计，无基因过滤。（b）基于GLM的似然比统计和t吨-统计数据显示，在去除大脑或UHR中读取少于20个的基因之前和之后。在这两个图中，如果一个基因的qRT-PCR绝对对数比小于0.2，则该基因被宣布为非DE基因，如果其绝对对数比大于2.0，则被宣布为DE基因。注意，根据mRNA-Seq和qRT-PCR，我们要求真阳性在同一方向上差异表达（见表1和方法）。

由于不同的mRNA-Seq DE测试显示出类似的行为，我们将从这里开始只关注基于GLM的似然比测试的结果。当使用不同的测试统计数据时，结果不会改变，除非已经注意到t吨-低计数基因的统计。

技术效果对差异表达的影响

差异表达中的基因长度偏差

根据mRNA-Seq分析，较长的转录物比相同水平表达的较短转录物产生更多“可排序”片段。基因计数和长度之间显然存在正相关，这种关联并没有通过基因长度的缩放完全消除，如在RPKM中[7]（[附加文件1：补充图S4]）。这表明较长基因的表达较高，或者基因计数对长度的非线性依赖。

如所述[11]，基因计数对长度的依赖性在信使核糖核酸序列DE结果中产生了“基因长度相关的偏差”：较长的基因往往具有更显著的DE统计数据（图2). 这里评估的所有mRNA-Seq DE统计数据都与读取计数的估计标准误差有内在的相关性。这在创建差异表达的“基因列表”方面是一个严重的缺陷，因为与具有较大影响的短基因相比，生成的列表可能倾向于具有较小潜在影响的长基因。仅考虑估计的折叠变化是不够的，因为这忽略了给定折叠变化和基因长度的相当大的标准误差范围。

图2

差异表达式统计，按长度.DE统计等级与UI基因长度的箱线图，UI基因长度至少为250 bp-long，在大脑和UHR中均为非零计数。（a） 德尔塔法t吨-统计数据。（b） Delta方法t吨-由基因长度平方根的倒数加权的统计学。

人们可以使用每个基因的固定碱基数来纠正DE统计的长度依赖性；通过从每个基因中随机选择250 bp重复DE分析，消除了DE显著性和长度之间的关联（[附加文件1：补充图S5]）。这也表明关联的原因是基因的长度，而不是较长基因的潜在表达水平的差异。然而，定长分析并不令人满意，因为它丢弃了大量数据，并且没有自然选择通用长度。

基于基因长度的加权分析可能构成对长度无关的DE滤波器的合理折衷。确实，缩放每个t吨-长度平方根的倒数统计提供了长度相关性排名（图2). 然而，选择截止点的问题仍然存在。根据下列假设[11]，与未加权t吨-统计数据和使用相同的跨基因截止值，对于给定的DE水平，功率随着基因长度的增加而增加。在相同的场景下，对于加权t吨-统计数据表明，I型错误率和功耗都随着长度的增加而降低。

底座校准方法的影响

在每个流细胞中保留八分之一的通道来对噬菌体phi X基因组DNA进行测序的做法在样本大小和平衡方面对实验设计具有重要意义。我们发现，在三个映射严格程度的每一个水平上，使用自动校准比使用标准phi X校准映射到基因组的读取数更多（图三). 纯度过滤完全匹配（FPM）读取不太可能包含排序错误，可以作为完全准确读取的代理。类似地，具有0、1或2不匹配（FMM）的纯过滤读取由FPM读取以及表示序列错误的读取组成。然后，可以将比率（FMM-FPM）/FMM视为测序错误率的粗略估计，假设没有SNP。对于所有车道，自动校准方法产生的错误率略低（约5%）。

图3

基线校准方法对读数映射的影响.四种生物样品（脑、UHR）和流-细胞（F2、F3）组合中每一种的每一条通道的平均读取计数条形图，包括φX校准和不带φX校准。读取分为三个嵌套类别：纯净过滤完全匹配读取（FPM）；具有0、1或2个不匹配的纯过滤读取（FMM）；0、1或2不匹配的未筛选读取（MM）。

阅读次数的增加在整个转录组中分布不均。大多数UI基因在校准方法之间的读数计数没有变化，而当使用自动校准时，大约25%的基因具有4个或更多的额外读数。当计算phi X和自动校准的每个基因的（FMM-FPM）/FMM比率时，自动校准平均产生约3.8%的错误率。

通过将观察到的差异与随机交换14条车道中每个车道的自动校准和φX校准读取计数集获得的差异的排列分布进行比较，评估两种校准方法之间表达测量差异的显著性。我们发现，在绝对表达测量方面，两种校准方法之间存在微小但具有统计意义的差异。然而，DE分析中使用的相对表达式度量似乎没有显著差异（参见[附加文件2：补充章节S8]）。

尽管我们的评估仅基于两个流动单元，但很明显，自动校准是有利的，因为它可以产生更平衡的设计，每个流动单元释放一个通道，并且每个通道产生更多更高质量的读数。

车道、流细胞和库准备效果

泊松分布已被证明在通过总车道数进行归一化后，能够很好地拟合复制车道上的基因级计数分布[4,6]；我们在MAQC数据和未发表的数据集方面的经验黑腹果蝇支持这一结论。泊松模型在不同生物体和不同测序设施中的良好性能有力地支持了其作为车道变化模型的有效性，并证明了通过求和跨车道合并读取计数的合理性。然而，请注意，在分析时，泊松分布的适用性值得怀疑生物复制（即，来自给定生物组中不同个体的样本，例如，患有相同类型癌症的患者）。SAGE文献中描述的负二项或经验贝叶斯方法的使用[21,22]，在这种可变性增加的环境中可能是合理的。

我们的分析也证实了前面提到的流式细胞之间的微小技术差异[6]尽管有证据表明流细胞之间的差异略大于复制通道之间的差异（[附加文件1：补充图S6c]）。无论其统计意义如何，估计的流细胞效应的大小都很小，因此只有在检测极小的生物效应时才会产生较小的影响；对于阅读次数超过3次的基因，几乎没有。

据我们所知，目前还没有对图书馆准备期间引入的技术变化进行公开审查；复制库准备工作既昂贵又耗时。库准备对总读取次数有明显影响（[附加文件1：补充图S1]）。在对总车道数的差异进行调整后，有证据表明，与流细胞和车道相比，复制库准备之间的差异增加了（[附加文件1：补充图S6d]）；然而，这种增加的变异性主要是由于高计数基因具有很强的检测小差异的能力。由于实验设计的原因，无法将文库制备效果与流动细胞和生物效果进行直接比较，但对估计差异大小的比较表明，文库制程效果远小于Brain和UHR之间的生物效果（图4)对于某些基因来说，比流细胞效应稍大（图4和[附加文件1：补充图S6]）。

图4

生物、文库制备和流细胞效应的比较UHR与脑生物效应的估计对数变化箱线图（[附加文件]中的GLM 22：补充表S4]）、针对生物学的流细胞效应调整（GLM 4）、流细胞内的库制备效应（GLM 7）。给出了总计数（黑色）和上四分位（蓝色）标准化的估计值。

大脑和UHR样本之间的生物学差异可能比通常观察到的要大得多；因此，变化的技术来源不一定总是无关紧要的。最后，我们注意到MAQC数据有点“理想”，即：（1）商业颗粒RNA已测序，（2）测序由Illumina公司内部执行。典型的信使核糖核酸序列实验始于从生物样本中提取RNA，提取过程中诱导的变异性可能比这里看到的技术变异性大得多。

mRNA-Seq数据的规范化

由于通道之间的读取总数不同，因此必须规范化读取计数，以允许跨通道或样本比较表达式度量。虽然这一主题在mRNA-Seq文献中受到的关注相对较少，但通常的做法是通过车道总数来衡量基因计数[6,7]. 然而，我们发现，更通用的基于分位数的程序与qRT-PCR产生更好的一致性，并且有望比单个家政基因的标准化更稳健。

在这里，我们评估了各种标准化程序，并重点关注两个主要问题：（1）标准化是否提高了DE检测（灵敏度）？（2） 标准化是否会导致重复之间的低技术可变性（特异性）？为了评估DE检测，我们依赖以下的qRT-PCR数据[13]作为判断真假阳性的金标准。由于qRT-PCR数据中的非DE基因数量有限，我们还评估了复制车道的泊松模型（附加文件中的GLM 12：补充表S4]）。

最简单的标准化形式是通过缩放车道中的基因计数来实现的我，通过单个车道特定因素d日 _我 本质上，这些全球的比例因子定义了无差异表达的无效假设：如果基因在车道上的计数比例与向量确定的比例相同d日 _我 ’s，那么它被认为是非差异表达的。

如前所述，标准总计数标准化可降低车道、流细胞和库准备之间的差异。以前还不清楚的是，这种标准化技术反映了相对较少的高计数基因的行为：5%的基因约占大脑和UHR总计数的50%。这些基因不能保证在不同的生物条件下具有相似的表达水平，并且在MAQC-2数据集的情况下，与大多数基因相比，它们在大脑中明显过度表达（图5).

图5

高表达基因的影响（a）脑（绿色）和UHR（紫色）样本总读取计数的累计百分比，从带有最高读取计数（跨越七个Brain或UHR通道）。累计读取计数标记为5%、10%、20%和30%的高表达基因。（b） 非正常计数的UHR/Brain表达fold-change的运行值，从带有最低的所有14条车道的总数。水平线对应于：所有基因的计数比率（黑色）、POLR2A基因的计数比例（红色）以及在至少一条通道中读取的基因的计数百分比上限四分位的比率（蓝色）。

因此，总计数归一化的性能对于检测DE来说并不是特别令人印象深刻（图6)：与微阵列数据相比，敏感性仅稍高，即使是表达差异相对较大的基因（>2绝对对数比）。当包含差异表达水平较低的基因（绝对对数比>0.5）时，性能并不比微阵列好（也许稍差）。这与一般的预期相矛盾，因为信使核糖核酸序列数据的噪声较小，因此更善于检测小的表达差异。对于较小水平的DE，使用总计数归一化估计的对数比率中的偏差使得测序估计不太准确。

我们评估了mRNA-Seq数据归一化的两种备选方案。一种方法依赖于一个单一的看家基因，如POLR2A，这是一种标准化qRT-PCR表达测量的技术。然而，这通常不是一个可行的解决方案，因为它是未知的先验的哪些基因具有稳定的表达水平（in[13]，POLR2A是在对多个平板上的UHR和Brain进行多次重复检测后选择的）。

图6

mRNA-Seq和微阵列差异表达调用qRT-PCR的比较：ROC曲线.所有三个平台通用的基因，用于qRT-PCR和测序（参见[附加文件2：补充章节S6]）。如果qRT-PCR绝对对数比（a）大于2或（b）大于0.5，则评估并申报DE；如果基因的绝对对数比小于0.2，则该基因被宣布为非DE基因。测序数据采用基于GLM的似然比检验。为mRNA-Seq提供了两种归一化程序：总计数（黑色）和上四分位（蓝色）归一化。使用RMA（灰色）对微阵列数据进行标准化。注意，根据mRNA-Seq和qRT-PCR，我们要求真阳性在同一方向上差异表达（见表1和方法）。

与标准的微阵列数据归一化技术类似，我们建议根据分位数等参数匹配基因计数的面间分布。例如，可以简单地按中间值缩放车道内的计数。在我们的案例中，由于零和低计数基因的优势，中位数对不同水平的测序工作没有信息。相反，我们在排除了所有通道中零读数的基因后，使用了per-lane上四分位（第75个百分位）（参见方法）。

与总计数标准化相比，POLR2A和上四分位数标准化显著降低了DE相对于qRT-PCR的偏差（图7和[附加文件1：补充图S7]），上四分位数偏差接近零。ROC曲线表明，在检测DE时，上四分位数和POLR2A归一化都明显优于总计数归一化（图6和[附加文件1：补充图S8a]），并导致测序相对于微阵列数据的灵敏度提高（图6和[附加文件1：补充图S9]）。

图7

mRNA-Seq和微阵列差异表达测量与qRT-PCR的比较差异散点图比较了qRT-PCR与（a）mRNA-Seq（使用标准总计数归一化）和（b）微阵列（使用标准RMA归一化。所示为所有三个平台之间共享的基因，根据mRNA-Seq和qRT-PCR都存在于大脑和UHR中（参见[附加文件2：补充部分S1]），且绝对qRT-PCR表达对数比小于4。（a）中的水平线表示测序数据经过标准总计数归一化（黑色）、POLR2A归一化（红色）、分位数归一化处理（黄色）、上四分位数归一处理（蓝色）后的中位数UHR/Brain log-ratio；（b）中的水平线显示了标准RMA标准化（黑色）和POLR2A标准化（红色）后微阵列数据的中位数UHR/Brain log-ratio。

仔细研究不同归一化程序的技术差异表明，与总计数归一化相比，上四分位归一化并没有显著增加可变性水平；POLR2A标准化的变化稍大，但仍具有可比性（图8).

最后，像在微阵列实验中经常做的那样，跨通道执行分位数归一化也是可行的[23]. 然而，这一更复杂的正常化战略似乎没有额外的好处。分位数归一化在ROC分析中的表现类似（图[附加文件1：补充图S8a]），并与上四分位归一化（图8). 我们再次回忆起MAQC数据的某种人为性质，这些数据基本上是在同一时间由一个实验室使用理想的RNA样本获得的。随着越来越多的数据可用，基因计数分布可能会有更大的变化，需要更积极的标准化。

图8

归一化过程的比较：泊松模型的拟合优度.乘法泊松模型（[附加文件]中的GLM 12：补充表S4]）在（a）总数、（b）POLR2A、（c）上四分位数和（d）分位数归一化后适用于MAQC-2实验中的七个脑区。计算并显示了优良率统计数据χ²分位数-分位数图。在基因的最高分位数中具有良好的fit统计信息的基因χ²-使用彩色绘图符号显示分布：红色（1,5]%、紫色（.1,1]%、金色[0,1]%。UHR的类似绘图显示相同的模式。

结论

我们的主要新发现是归一化对差异表达结果的影响程度：归一化过程之间的敏感性差异比测试统计数据之间的差异更大。尽管标准总计数归一化会导致复制车道之间的泊松变化，但当以qRT-PCR为基准时，它的检测灵敏度较差。相反，我们建议通过基因计数分布的分位数（上四分位数）来缩放基因计数，并表明这种标准化可以提高灵敏度，而不会丢失特异性。

MAQC数据集的一个重要方面是大脑和UHR之间基因表达的巨大差异，这可能会影响对所述分析的解释。通常，基因表达分析考虑更密切相关的样本集，预计只有相对较少的基因会有差异表达。相比之下，在Brain和UHR的比较中，通过qRT-PCR检测的基因只有5-30%被认为是不-差异表达（取决于用于校正第页-值）。事实上，qRT-PCR实验可能没有真正的非DE基因，而是每个基因的表达差异非常小。这在指定一组真负值时可能会产生错误；我们试图通过仔细严格地定义真阴性，并通过评估定义变化的影响来控制这种情况（参见[附加文件2：补充章节S6]）。

此外，大脑和UHR样本之间转录谱的极端差异意味着第页-如果所有基因都是非DE基因，测序实验得出的值会小于预期值第页-非DE基因的值（根据qRT-PCR）不遵循预期的均匀分布，但明显向零偏移（[附加文件1：补充图S10]）。微阵列数据显示了相同的行为（[附加文件1：补充图S10]），表明这是由考虑中的样本引起的，而不是由统计方法的固有问题引起的。与差异表达的qRT-PCR测试相反，用于测序数据的测试考虑了每个基因的总读取数，因此，对于较长的基因往往具有更大的能力。

另一个可能的批评是，UQ相对于总计数标准化的改进是由于这种标准化与归一化使用qRT-PCR数据而不是正确反映实际生物差异的程序；换言之，UQ标准化可能与POLR2A除法的效果密切匹配，就像qRT-PCR数据所做的那样，但不是潜在的生物学。事实上，与标准微阵列分位数归一化相比，POLR2A对微阵列数据的额外缩放略微改善了ROC（[附加文件1：补充图S8b]）。然而，更可能的是，由于总计数正常化依赖于高计数基因，因此不能很好地反映生物差异。通过仔细观察POLR2A基因可以看出这一点，该基因被选为qRT-PCR数据的参考，因为其在UHR和大脑中的表达在许多qRT-PCR重复物中非常相似[13]：对于总计数归一化，估计POLR2A的UHR到脑折叠变化为1.3，而对于四分位上归一化为0.97，对于微阵列数据为0.90。

关于DE检验统计，基于GLM的似然比统计和Fisher精确统计在低计数基因的敏感性和处理方面表现同样好。我们发现似然比测试因其通用性而具有吸引力。实际上，使用GLM框架，可以调整潜在的混淆变量，包括定量协变量，例如样本年龄，以及适应不同的计数分布（过度分散情况下为负二项式）。

这里考虑的所有DE方法的一个严重问题是读计数对功率的固有依赖性，而读计数又与基因表达水平和长度有关。由于大多数DE研究产生的基因列表通常与功能注释（例如GO）相关，因此不希望由长度等特征驱动显著性值。基于基因长度的加权分析可能会导致基因的合理长度依赖性排序，这将使具有较大影响的短基因比具有较小影响的长基因显著增加。

我们发现，车道和流动中心之间的技术差异很小，而图书馆准备工作之间的技术变化稍大。然而，在所有情况下，对差异表达结果的影响都很小。如上所述，MAQC数据集是不寻常的，因为我们预计Brain和UHR之间的表达差异会非常大，并且由于RNA的高质量，仅会产生很小的库准备效应。实际上，库准备效应可能更接近于生物效应。

我们已经证明，虽然在分析管道的早期阶段，φX和自动校准之间存在一些差异，但差异表达方面的差异很小。总的来说，自动校准似乎是有利的，因为它产生了更平衡的设计，每个流单元释放一条通道，并且每个通道产生更多更高质量的读取。

这项工作以及其他工作中进行的分析都是基于表达谱的“全基因”观点。我们使用非常受限的UI基因定义评估了技术效果、φX校准和标准化方法。我们将自己局限于这样一个严格的定义，以确保评估不受选择性剪接或重叠基因的影响。我们的UI基因定义过于简单，因为大量的生物学相关信息丢失了；我们排除了50%以上属于集合基因的读取。

随着高通量测序变得越来越普遍，我们精确描述样本转录组的能力将显著增强。更精细的分析，如异构体水平的表达、等位基因特异性表达和基因组注释（分割），涉及比较样本内不同区域与样本间相同区域的差异。此类分析需要了解序列组成对基础覆盖率的影响，以解释基础水平计数分布的异质性