Differential expression analysis for sequence count data

Anders, Simon; Huber, Wolfgang

doi:10.1186/gb-2010-11-10-r106

方法
开放式访问
出版：2010年10月27日

序列计数数据的差异表达分析

基因组生物学 体积 11，物品编号：106兰特(2010)引用这篇文章

38.1万访问
10公里引文
107海拔高度
韵律学细节

摘要

高通量测序分析，如RNA-Seq、ChIP-Seq或条形码计数，以计数数据的形式提供定量读数。为了正确推断此类数据中的差分信号并具有良好的统计能力，需要估计整个动态范围内的数据可变性和合适的误差模型。我们提出了一种基于负二项分布、方差和均值通过局部回归联系的方法，DESeq公司，作为R/生物导体包。

背景

DNA片段的高通量测序用于一系列定量分析。这些分析的一个共同特点是，它们对大量DNA片段进行测序，这些DNA片段反映了例如生物系统的RNA分子库（RNA-Seq[1,2])或核苷酸结合分子的DNA或RNA相互作用区域（ChIP-Seq[三]、HITS-CLIP[4]). 通常，这些读取是根据它们到目标基因组公共区域的映射分配给一个类别的，在RNA-Seq的情况下，每个类别代表一个目标转录物，在ChIP-Seq的情况中，每个类别表示一个结合区域。一个重要的汇总统计是一个类中的读取次数；对于RNA-Seq，这个读取计数已经发现（很好地近似）与目标转录物的丰度线性相关[2]. 兴趣在于比较不同生物条件下的读取计数。在最简单的情况下，逐个类单独进行比较。我们将使用该术语基因与类同义，即使类也可能指转录因子结合位点，甚至条形码[5].

我们希望使用统计测试来确定，对于给定的基因，读取计数中观察到的差异是否显著，也就是说，是否大于仅仅由于自然随机变化而预期的值。

如果读取数据是从具有给定的固定基因片段的人群中独立采样的，则读取计数将遵循多项式分布，该分布可近似为泊松分布。

因此，泊松分布被用来测试微分表达式[6,7]. 泊松分布有一个单一的参数，该参数由其平均值唯一确定；它的方差和所有其他属性都是由此而来的；特别是，方差等于平均值。然而，已经注意到了[1,8]泊松分布的假设过于严格：它预测的变化比数据中看到的变化要小。因此，由此产生的统计测试并不能像广告中那样控制I类错误（错误发现的概率）。我们将在稍后的讨论中展示这方面的实例。

为了解决这个所谓的过度分散问题，有人建议用负二项（NB）分布对计数数据建模[9]，此方法用于边缘RSAGE和RNA-Seq分析包[8,10]. NB分布具有由平均值唯一确定的参数μ和方差σ²然而，感兴趣的数据集中的重复次数往往太少，无法可靠地估计每个基因的平均值和方差参数。对于边缘RRobinson和Smyth假设[11]平均值和方差之间的关系σ²=μ+αμ²，具有单个比例常数α这在整个实验中都是一样的，并且可以从数据中进行估计。因此，每个基因只需估计一个参数，就可以应用于少量重复的实验。

在本文中，我们通过允许更一般的、数据驱动的方差和均值关系来扩展此模型，提供了一种有效的算法来将模型拟合到数据，并表明它提供了更好的拟合（第模型). 因此，可以在数据的动态范围内更均衡地选择差异表达的基因（第差分表达式测试). 我们通过将其应用于四个数据集来演示该方法（第节应用)并讨论它与其他方法的比较（第节结论)。

结果和讨论

模型

描述

我们假设样本中的读取次数j个分配给基因的我可以用负二项（NB）分布建模，

{K（K）}_{我 j个} ~ 注 (μ_{我 j个}, σ_{我 j个}^{2}),

（1）

它有两个参数，平均值μ_伊吉以及方差 $σ_{我 j个}^{2}$ .读取计数K（K）_伊吉是非负整数。分布概率在补充注释A中给出。（所有补充注释都在附加文件中1.）当存在过度分散时，NB分布通常用于建模计数数据[12].

实际上，我们不知道参数μ_伊吉和 $σ_{我 j个}^{2}$ ，我们需要根据数据进行估算。通常，重复次数很少，需要进行进一步的建模假设，以获得有用的估计。在本文中，我们开发了一种基于以下三个假设的方法。

首先，平均参数μ_伊吉即基因的观察计数的期望值我在样品中j个，是条件相关的per-gene值的乘积q个_我,_ρ(_j个)（其中ρ(j个)是样品的实验条件j个)和尺寸系数秒_j个,

μ_{我 j个} = {q个}_{我, ρ {(j个)}^{秒} j个} .

(2)

q个_{i、 ρ(j个)}与基因片段真实（但未知）浓度的期望值成正比我在条件下ρ(j个). 尺寸系数秒_j个表示库的覆盖范围或采样深度j个，我们将使用该术语普通比例尺对于数量，例如q个_我,ρ(j个)，通过除以进行覆盖调整秒_j个.

第二，方差 $σ_{我 j个}^{2}$ 是a的总和散粒噪声项和a原始方差项,

σ_{我 j个}^{2} = \underset{射击 噪音}{\underset{⏟}{μ_{我 j个}}} + \underset{未经加工的 方差}{\underset{⏟}{秒_{j个}^{2} {v（v）}_{我, ρ (j个)}}} .

(3)

第三，我们假设per-gene原始方差参数v（v）_我,ρ是的平滑函数q个_我,ρ,

{v（v）}_{我, ρ (j个)} = {v（v）}_{ρ} ({q个}_{我, ρ (j个)}) .

(4)

这种假设是必要的，因为复制次数通常太少，无法精确估计基因的方差我根据这个基因的现有数据。这种假设使我们能够汇集来自表达强度相似的基因的数据，以进行方差估计。

方程中方差的分解(三)由以下层次模型驱动：我们假设基因片段的实际浓度我在样品中j个与随机变量成比例R（右）_伊吉这样基因片段的比率我已排序为秒_j个第页_伊吉对于每个基因我和所有样品j个条件的ρ，的R（右）_伊吉i.i.d.是平均值q个_iρ和方差v（v）_iρ因此，计数值K（K）_伊吉，条件为R（右）_伊吉=第页_伊吉，是泊松分布率秒_j个第页_伊吉.的边际分布K（K）_伊吉-当允许变化时R（右）_伊吉-具有平均值μ_伊吉和（根据总方差定律）方程式中给出的方差(三). 此外，如果R（右）_伊吉根据伽马分布建模K（K）_伊吉为NB（例如，参见[12]，第4.2.2节）。

适合的

我们现在描述如何将模型拟合到数据中。这些数据是n个×米计数表，k个_伊吉，其中我= 1,...,n个索引基因，以及j个= 1,...,米为样本编制索引。该模型有三组参数：

（i）米尺寸系数秒_j个; 样本中所有计数的期望值j个与…成比例秒_j个.

（ii）对于每个实验条件ρ,n个表达式强度参数q个_iρ; 它们反映了基因片段的预期丰度我处于条件下ρ也就是基因计数的期望值我与…成比例q个_iρ.

（iii）平滑功能v（v）_ρ:ℝ⁺→ℝ⁺; 对于每个条件ρ,v（v）_ρ对原始方差的相关性建模v（v）_iρ关于期望均值q个_iρ.

尺寸系数的目的秒_j个是从不同样本中渲染计数，这些样本可能已按不同深度排序，具有可比性。因此，比率( $E类$ K（K）_伊吉)/( $E类$ K（K）_ij’)相同基因的预期计数我在不同的样品中j个和j’应等于尺寸比秒_j个/秒_j’if基因我没有差异表达或样本j个和j’是重复的。读取总数，∑_我 k个_伊吉，似乎是测序深度的一个很好的衡量标准，因此是秒_j个然而，实际数据的经验表明，情况并非总是如此，因为一些高度差异表达的基因可能会对总读取计数产生强烈影响，导致总读取计数的比率不能很好地估计预期计数的比率。

因此，为了估计规模因素，我们取观察计数比率的中位数。将刚刚概述的程序推广到两个以上样本的情况，我们使用：

{\hat{秒}}_{j个} = \underset{我}{中值的} \frac{{k个}_{我 j个}}{{(\prod_{v（v） = 1}^{米} {k个}_{我 v（v）})}^{1 / 米}} .

(5)

该表达式的分母可以解释为通过对样本取几何平均值获得的伪参考样本。因此，每个尺寸系数估计值 ${\hat{秒}}_{j个}$ 计算为j个-第个样本的计数与伪参考的计数相同。（注：在审查这份手稿时，罗宾逊和奥什拉克[13]建议了类似的方法。）

估计q个_iρ，我们使用样本计数的平均值j个与条件相对应ρ，转换为普通规模：

{\hat{q个}}_{我 ρ} = \frac{1}{米_{ρ}} \sum_{j个 : ρ (j个) = ρ} \frac{{k个}_{我 j个}}{{\hat{秒}}_{j个}},

(6)

哪里米_ρ是条件的重复次数ρ和在这些重复上运行。功能v（v）_ρ，我们首先计算普通尺度上的样本方差

{w个}_{我 ρ} = \frac{1}{米_{ρ} - 1} \sum_{j个 : ρ (j个) = ρ} {(\frac{{k个}_{我 j个}}{{\hat{秒}}_{j个}} - {\hat{q个}}_{我 ρ})}^{2}

(7)

并定义

{z（z）}_{我 ρ} = \frac{{\hat{q个}}_{我 ρ}}{米_{ρ}} \sum_{j个 : ρ (j个) = ρ} \frac{1}{{\hat{秒}}_{j个}} .

(8)

在附加文件的补充注释B中1我们证明了这一点w个_iρ-z（z）_iρ是原始方差参数的无偏估计量v（v）_iρ方程式的(三)。

然而，对于少量重复，米_ρ与应用程序中的典型情况一样w个_iρ高度可变，并且w个_iρ-z（z）_iρ对于统计推断来说，这不是一个有用的方差估计量。相反，我们使用局部回归[14]关于图 $({\hat{q个}}_{我 ρ}, {w个}_{我 ρ})$ 获得平滑函数w个_ρ(q个)，使用

{\hat{v（v）}}_{ρ} ({\hat{q个}}_{我 ρ}) = {w个}_{ρ} ({\hat{q个}}_{我 ρ}) - {z（z）}_{我 ρ}

(9)

作为我们对原始方差的估计。

需要注意避免局部回归中的估计偏差。w个_iρ是平方随机变量和残差的总和 ${w个}_{我 ρ} - w个 ({\hat{q个}}_{我 ρ})$ 倾斜。以下参考文献[15]，第8章和[14]，第9.1.2节，我们使用伽马族的广义线性模型进行局部回归，使用位置匹配包装[16].

差异表达测试

假设我们有米_A类复制生物条件A和米_B类条件B的样本。每个基因我，我们想权衡数据中该基因在两种条件下差异表达的证据。特别是，我们想测试虚假设q个_国际机场=q个_国际银行，其中q个_国际机场是条件A和q的样本的表达式强度参数_国际银行对于条件B，为此，我们将每个条件中的总计数定义为测试统计，

\begin{matrix} {K（K）}_{我 A类} = \sum_{j个 : ρ (j个) = A类} {K（K）}_{我 j个}, & {K（K）}_{我 B类} = \sum_{j个 : ρ (j个) = B类} {K（K）}_{我 j个}, \end{matrix}

(10)

和他们的总金额K（K）_是=K（K）_国际机场+K（K）_国际银行根据上一节中描述的错误模型，我们在下面显示，在零假设下，我们可以计算事件的概率K（K）_国际机场=一和K（K）_iB公司=b条对于任何一对数字一和b条。我们将此概率表示为第页(一,b条). 这个P（P）一对观测计数和的值(k个_国际机场,k个_国际银行)那么所有概率的总和小于或等于第页(k个_国际机场,k个_国际银行)，假设总金额为k个_是:

{第页}_{我} = \frac{\sum_{\begin{matrix} 一 + b条 = {k个}_{我 秒} \\ 第页 (一, b条) \leq 第页 ({k个}_{我 A类} {k个}_{我 B类}) \end{matrix}} 第页 (一, b条)}{\sum_{一 + b条 = {k个}_{我 秒}} 第页 (一, b条)} .

(11)

变量一和b条在上述总和中，取0，。。。，k个_我秒到目前为止，提出的方法遵循了罗宾逊和斯迈思的方法[11]与其他条件测试类似，如Fisher精确测试。（参见参考[17]，第3章，讨论测试中条件反射的优点。）

计算第页(一,b条). 首先，假设在零假设下，不同样本的计数是独立的。然后，第页(一,b条)=优先级(K（K）_我A类=一)优先级(K（K）_我B类=b条). 因此，问题是计算事件的概率K（K）_我A类=一类似地，和K（K）_我B类=b条.随机变量K（K）_我A类是以下各项的总和米_A类

NB-分布随机变量。我们通过NB分布来近似其分布，NB分布的参数我们从K（K）_伊吉为此，我们首先根据两种条件的计数计算合并平均值估计值，

{\hat{q个}}_{我 0} = \sum_{j个 : ρ (j个) \in {A类, B类}} {k个}_{我 j个} / 秒_{j个},

(12)

这说明了虚假设规定q个_我A类=q个_我B类条件A的平均值和方差之和为

{\hat{μ}}_{我 A类} = \sum_{j个 \in A类} 秒_{j个} {\hat{q个}}_{我 0},

(13)

{\hat{σ}}_{我 A类}^{2} = \sum_{j个 \in A类} {\hat{秒}}_{j个} {\hat{q个}}_{我 0} + {\hat{秒}}_{j个}^{2} {\hat{v（v）}}_{A类} ({\hat{q个}}_{我 0}) .

（14）

附加文件中的补充注释C1描述了NB的分布参数K（K）_我A类可以根据 ${\hat{μ}}_{我 A类}$ 和 ${\hat{σ}}_{我 A类}^{2}$ （为了避免偏差，我们不直接匹配力矩，而是匹配不同的一对分布统计数据。）K（K）_我B类通过类比获得。

附加文件中的补充注释D1解释了我们如何计算方程式中的总和(11)。

应用

数据集

我们根据以下数据集给出了结果：

飞行胚胎中的RNA-Seq

B.Wilczynski、Y.-H.Liu、N.Delhomme和E.Furlong在苍蝇胚胎中进行了RNA-Seq实验，并在发表之前与我们分享了部分数据。在这个数据集的每个样本中，都有一个基因被设计成过表达，我们比较了两个生物复制品，每一个都是这样的条件，以下分别表示为“a”和“B”。

神经干细胞的标记序列

Engström公司等。[18]执行Tag-Seq[19]用于神经细胞的组织培养，包括四个来自胶质母细胞瘤衍生神经干细胞（GNS）和两个来自非癌神经干细胞。由于每种组织培养物来源于不同的受试者，因此具有不同的基因型，这些数据显示出高度的可变性。

酵母的RNA-Seq

那加拉克什米等。[1]对酿酒酵母文化。他们测试了两种库准备协议，数据传输和右侧，并对每个方案进行了三次测序，因此对于每个方案的第一次运行，他们有一个进一步的技术复制（相同的培养物，复制的文库制备）和一个进一步生物复制（不同的培养物）。

HapMap样本的ChIP-Seq

卡索夫斯基等。[20]通过ChIP-Seq比较了10名人类个体DNA区域的蛋白质占位。他们编制了聚合酶II和NF-κB的区域列表，并计算了每个样本对应于每个区域的读取数。这项研究的目的是调查不同地区的职业在个人之间的差异。

方差估计

我们首先演示方差估计。图1a个显示了样本方差w个_iρ（方程式(7))根据平均值绘制 ${\hat{q个}}_{我 ρ}$ （方程式(6))对于条件A类在飞行RNA-Seq数据中。图中还显示了局部回归拟合w个_ρ(q个)和散粒噪音 ${\hat{秒}}_{j个} {\hat{q个}}_{我 ρ}$ .在图中1亿，我们绘制了平方变异系数（SCV），即方差与均方的比率。在该图中，橙色线和紫色线之间的距离是由于生物采样引起的噪声的SCV（参见等式(三)).

图中的许多数据点1亿这远远高于拟合的橙色曲线，可能会导致局部回归拟合效果不佳。然而，预计残差分布会出现强烈的偏差。参见附加文件中的补充注释E1以了解适用于验证匹配性的诊断的详细信息和讨论。

测试

为了验证DESeq公司保持对I型错误的控制，我们对比了其中一个复制品的情况A类在对另一个样本的飞行数据中，对两个样本使用从两个重复中估计的方差函数。图2显示了P（P）通过比较获得的值。为了控制I类错误P（P）值低于阈值α必须≤α也就是说，ECDF曲线（蓝线）不应超过对角线（灰线）。如图所示，I类错误由边缘R和DESeq公司，但不是基于泊松的χ²测试。后者低估了数据的可变性，因此会导致许多假阳性拒绝。除了对真实数据进行评估外，我们还验证了DESeq公司对上述错误模型生成的模拟数据进行I类错误控制；参见附加文件中的补充注释G1接下来，我们对比了两者A类两个样本的对比B类样本。使用上一节中描述的过程，我们计算了P（P）每个基因的值。图三显示获得的褶皱变化和P（P）值。12%的P值低于5%。使用Benjamini和Hochberg程序进行多次试验的调整[21]864个基因（17605个）的错误发现率（FDR）为10%，差异表达显著。图中用红色标记三演示了检测差异表达的能力如何依赖于总计数。具体来说，低计数的强散粒噪声导致测试过程只调用非常高的折叠变化。还可以看出，对于大约100以下的计数，即使计数水平略有增加，也会降低散粒噪声的影响，从而降低折叠式变换的要求，而在较高计数时，当散粒噪声变得不重要时（参见图1亿)，折叠式截止线仅弱依赖于计数水平。这些图有助于指导实验设计：对于弱表达基因，在散粒噪声很重要的区域，可以通过更深入的测序来增加功率，而对于更高计数状态，只有通过进一步的生物复制才能实现功率的增加。

与edgeR的比较

我们还分析了数据边缘R（版本1.6.0[8,10,11]). 我们跑了边缘R有四种不同的设置，即在通用分散模式和分段分散模式下，或者使用根据DESeq公司或获取顺序读取的总数。结果与这些选择没有太大关系，这里我们报告了标签分散模式的结果DESeq公司-估计的规模因素。（复制本文中报告的所有分析、数字和数字所需的R代码在附加文件中提供2; 此外，本补充提供了以下其他设置的结果边缘R。原始数据可以在附加文件中找到三）

回到图1我们看到了边缘R的方差的单值离散估计值低于DESeq公司弱表达基因和强表达基因的表达量较高。因此，如图所示2边缘R对于低表达基因是反保守的。然而，它通过对强表达基因更加保守来弥补这一不足，因此，平均而言，I型错误控制得以维持。

然而，在不同条件之间的测试中，这种行为可能会导致发现列表中出现偏差；对于当前数据，如图4显示，弱表达基因似乎过表达，而很少有平均水平高的基因被称为差异表达边缘R虽然两种方法的总体敏感性似乎相当(DESeq公司报告了864次点击，边缘R1127次点击），DESeq公司产生的结果在动态范围内更加平衡。

神经干细胞数据也获得了类似的结果，该数据集具有不同的生物背景和不同的噪声特性（参见附加文件中的补充注释F1). 本研究中提出的方差估计方案的灵活性似乎在一系列应用中比现有方法提供了真正的优势。

在没有复制的情况下工作

DESeq公司允许在一种或两种条件下分析没有生物复制的实验。虽然人们可能不想从这样的分析中得出强有力的结论，但它可能仍然对探索和假设生成有用。

如果复制仅适用于其中一种条件，则可以选择假设根据该条件的数据估计的方差-均值相关性也适用于未复制的条件。

如果这两种情况都没有重复，我们仍然可以基于以下假设进行分析：对于大多数基因，没有真正的差异表达，并且可以通过将两个样本视为重复来估计有效的均值-方差关系。少数差异丰富的基因将充当异常值；然而，它们不会对gamma-family GLM拟合产生严重影响，因为形状参数低值的gamma分布具有很重的右尾。预计可能会高估方差，这将使该方法变得保守。

我们使用飞行RNA-Seq和神经细胞Tag-Seq数据进行了这样的分析，将这两个数据集限制为两个样本，每个条件一个样本。对于神经细胞数据，正如预期的那样，估计的方差函数略高于从全球导航系统和NS公司复制。

使用它来测试差异表达，在FDR=10%时仍然发现269次点击，其中202次是所有可用样本的更可靠分析的612次点击中的一次。然而，就飞行RNA-Seq数据而言，862次命中中只有90次（11%）被恢复（有两次新命中）。这些观察结果的解释是，在神经细胞数据中，重复之间的变异性并不比条件之间的变异小得多，这使得后者成为前者的可用替代品。另一方面，对于苍蝇数据，复制之间的变异性远小于条件之间的变异，这表明复制提供了有关数据中实验变异的重要信息，但在其他方面没有可用信息（另见下一节）。

方差稳定变换

给定方差-均值依赖性，方差稳定变换（VST）是一个单调映射，使得对于变换后的值，方差（近似）与均值无关。使用方差-均值相关性w个(q个)估算人DESeq公司，VST由下式给出

τ (κ) = {¦Β}^{κ} \frac{d日 q个}{\sqrt{w个 (q个)}} .

(15)

应用转换τ到通用刻度计数数据，k个_伊吉/秒_j个，生成在整个动态范围内方差大致相同的值。VST的一个应用是样本聚类，如图所示5; 这种方法比在未转换的计数数据上定义合适的距离度量更简单，其选择并不明显，并且可能不容易与可用的聚类或分类算法结合（这些算法往往是为具有类似分布特性的变量设计的）。

ChIP-Seq公司

DESeq公司也可用于分析比较ChIP-Seq分析。卡索夫斯基等。[20]分析HapMap个体的转录因子结合，并计算每个样本中映射到预定结合区域的读取次数。我们从他们的数据集中考虑了两个个体，即HapMap ID GM12878和GM12891，这两个个体至少进行了四次重复，并测试了这两个区域的不同职业。图的左上两个面板6显示同一个体内的比较，表明I型错误由以下因素控制DESeq公司。使用Benjamini-Hochberg调整，在10%FDR时无显著区域。然而，当对比这两个个体时，发现了不同的职业，在19028个区域中，有4460个在仅使用两个重复时显著，而在使用四个重复时则有8442个显著（右上方两个面板）。

使用另一种方法，卡索夫斯基等拟合泊松族的广义线性模型（GLM）。这（图的下一行6)导致了P（P）即使是在同一个体内进行比较时，也会出现值，这表明泊松GLM低估了方差，而P值的直接使用会导致反保守（过度乐观）偏差。卡索夫斯基等解决了这一问题，并通过使用其他标准来调用差异职业来调整偏差。

结论

为什么有必要为序列计数数据开发新的统计方法？如果有大量重复数据可用，则可以通过使用非参数方法，如基于等级或排列测试，来避免数据分布问题。然而，最好（也可能）考虑在每个条件下使用较少重复次数的实验。为了将观察到的差异与预期的随机变化进行比较，我们可以通过两种方式改进对后者的描述：首先，我们可以使用分布族，例如正态分布、泊松分布和负二项分布，以确定微分表达式统计的高阶矩，从而确定尾部行为，基于观察到的低阶矩，如均值和方差。其次，我们可以基于来自不同基因的数据遵循相似的变异模式的概念，在基因之间共享信息，例如分布参数。在这里，我们描述了这种方法的一个实例，现在我们将讨论我们所做的选择。

分销的选择

虽然对于大型计数，正态分布可能提供了复制间变异性的良好近似值，但对于较低的计数值，情况并非如此，其离散性和偏度意味着从正态近似值计算的概率估计值将不充分。

对于泊松近似，Marioni的工作是一篇关键论文等。[6]，研究了技术的RNA-Seq的再现性。他们从两个组织样本中提取了总RNA，一个来自同一个人的肝脏，一个来自同一个人的肾脏。他们从每个RNA样本中取出七份样本，根据Illumina推荐的方案从每个样本中制备一个文库，并在Solexa基因组分析仪的一条通道上对每个文库进行取样。对于每个基因，他们计算了来自同一组织样本的七个计数的方差，发现与泊松模型预测的方差非常吻合。与我们在第节中的论点一致模型泊松散粒噪声是计数过程中预期的最小变化量。因此，马里奥尼等得出结论，RNA-Seq的技术重复性很好，技术重复之间的差异接近散粒噪声极限。从这个有利位置来看，马里奥尼等.（和类似的Bullard等。[22])建议使用泊松模型（以及Fisher精确检验或似然比检验作为近似值）来测试基因在两个样本之间是否存在差异表达。值得注意的是，拒绝此类测试只会告诉我们，两个样本中平均计数之间的差异大于预期值技术的复制。因此，我们不知道这种差异是否是由于不同的组织类型（肾脏而非肝脏）造成的，也不知道如果比较同一肝脏不同部位的两个样本或两个人的肝脏，是否会发现同样大小的差异。

图1结果表明，散粒噪声仅在非常低的计数值中占主导地位，而对于中等计数，样本之间生物变化的影响超过散粒噪声几个数量级。

通过技术复制品与生物复制品的比较证实了这一点[1]. 在图中7我们用过DESeq公司获得Nagalakshmi数据的方差估计等。[1]. 分析表明，技术复制之间的差异几乎不超过散粒噪声水平，而生物复制之间的差别更大。基于泊松模型的微分表达式测试，如参考文献中讨论的那些[6,7,20,22,23]因此，应谨慎解释，因为它们可能严重低估了生物变异性的影响，尤其是对高表达基因的影响。

因此，最好使用允许过度分散的模型。对于泊松分布，方差和平均值相等，而负二项分布是一种泛化，允许方差更大。使用此发行版的最先进的已发布方法可能是边缘R[8].DESeq公司其基本思想归功于边缘R，但在几个方面有所不同。

基因间信息共享

首先，我们发现使用总读取计数作为测序深度的估计值，从而调整样本之间的观察计数（如Robinson建议的那样等。[8]和其他）可能会导致重复之间的明显差异，因此检测真实差异的能力较差。

DESeq公司使用更可靠的尺寸估算公式(5); 事实上，边缘R当提供这些尺寸估计值时，的功率会增加。（注：在审查本文件时，边缘R修改为使用Oshlack和Robinson的方法[13].)

对于实践中经常遇到的少量重复，不可能同时获得NB分布方差和平均参数的可靠估计。边缘R通过估计单个共同色散参数。在我们的方法中，我们利用这种可能性来估计更灵活、平均相关的局部回归。典型实验中可用的数据量足够大，足以对色散进行足够精确的局部估计。在RNA-Seq典型的大动态范围内，原始SCV经常出现显著变化，考虑到这一点，允许DESeq公司避免在差分表达式调用中偏向动态范围的某些区域（见图2和4)。

这种灵活性是DESeq公司和边缘R模拟结果表明边缘R和DESeq公司如果提供了具有恒定SCV的人工数据，则进行比较（附加文件中的补充注释G1)。边缘R试图通过允许用per-gene经验方差调整基于模型的方差估计来弥补单参数噪声模型的刚性。经验贝叶斯程序，类似于最初为利马包装[24–26]，决定了如何将这两种信息源最佳地组合在一起。然而，对于典型的低复制数，这种所谓的分段分散模式似乎没有什么影响（图4)甚至减少边缘R的功率（附加文件中的补充注释F1)。

第三，我们提出了一种简单而稳健的方法来估计数据的原始方差。罗宾逊和史密斯[11]采用了一种称为分位数调整条件最大似然的技术，以找到原始SCV的无偏估计。这个分位数调整指一种基于等级的过程，该过程修改数据，使数据似乎来源于相同库大小的样本。在DESeq公司，不同的库大小仅通过线性缩放来解决（方程(2)和(三))这表明分位数调整是一个不必要的复杂问题。我们为此付出的代价是，我们需要对方程中的NB变量之和进行近似(10)本身是NB分发的。虽然分位数调整和我们的近似值似乎都没有在实践中引起关注的理由，DESeq公司的方法计算速度更快，而且概念上可能更简单。

第四，我们的方法提供了有用的诊断。附加文件中的补充图S3等图2有助于判断测试的可靠性。在图中1亿和7很容易看出，在哪个平均值上，生物变异性比散粒噪声占主导地位；这些信息对于决定测序深度或生物复制次数是否是检测能力的限制因素很有价值，因此有助于规划实验。图中的热图5对数据质量控制很有用。

材料和方法

R包DESeq

我们将该方法作为统计环境R的包来实现[27]并在生物导体项目中分发[28]. 作为输入，它需要一个计数数据表。数据以及元数据，例如样本和基因注释，都是用S4类管理的计数数据集，源自电子设置，Bioconductor的标准数据类型，用于表格式数据。该包提供了执行分析的高级功能，如第节所示应用程序只需几个命令，就可以让对R知之甚少的研究人员使用它。这可以通过软件包附带的文档（软件包vignette）中的示例来证明。此外，还为希望偏离标准工作流程的高级用户提供了较低级别的功能。典型计算，如第节所示的分析应用，需要在个人计算机上花费几分钟的时间。

此处显示的所有分析都是用DESeq公司。希望对其进行详细检查的读者将找到一个Swave文档，其中包含分析代码的注释R代码作为附加文件2和附加文件中的原始数据三.

DESeq公司可从Bioconductor存储库中获得Bioconductor包[28]和来自[36].

缩写

ChIP-Seq公司：: 免疫沉淀染色质的高通量测序
ECDF公司：: 经验累积分布函数
财务总监：: 假发现率
GLM公司：: 广义线性模型
RNA-Seq：: RNA（高通量）测序
SCV公司：: 平方变异系数
注意：: 阴性肿瘤（分布）
VST（垂直变速箱）：: 方差稳定变换。

工具书类

Nagalakshmi U，Wang Z，Waern K，Shou C，Raha D，Gerstein M，Snyder M：通过RNA测序确定的酵母基因组转录图谱。科学。2008, 320: 1344–1349. 10.1126/科学1158441。
第条中国科学院谷歌学者
Mortazavi A、Williams BA、McCue K、Schaeffer L、Wold B：通过RNA-Seq对哺乳动物转录体进行定位和量化。自然方法。2008年，5:621–628。10.1038/nmeth.1226。
第条中国科学院谷歌学者
Robertson G、Hirst M、Bainbridge M、Bilenky M、Zhao Y、Zeng T、Eukilchen G、Bernier B、Varhol R、Delaney A、Thiessen N、Griffith OL、He A、Marra M、Snyder M、Jones S：使用染色质免疫沉淀和大规模平行测序的STAT1 DNA关联的基因组全谱。自然方法。2007, 4: 651–657. 10.1038/nmeth1068。
第条中国科学院谷歌学者
Licatalosi DD、Mele A、Fak JJ、Ule J、Kayikci M、Chi SW、Clark TA、Schweitzer AC、Blume JE、Wang X、Darnell JC、Darnill RB:HITS-CLIP产生了大脑替代RNA处理的全基因组见解。自然。2008年，456:464–469。10.1038/nature07488。
第条中国科学院谷歌学者
Smith AM、Heisler LE、Mellor J、Kaper F、Thompson MJ、Chee M、Roth FP、Giaever G、Nislow C：通过深度条形码测序进行定量表型分析。《基因组研究》，2009，19:1836–1842。10.1101/gr.093955.109。
第条中国科学院谷歌学者
Marioni JC、Mason CE、Mane SM、Stephens M、Gilad Y:RNA-seq：技术再现性评估和与基因表达阵列的比较。《基因组研究》2008，18:1509–1517。10.1101克/克079558.108。
第条中国科学院谷歌学者
Wang L，Feng Z，Wang X，Wang X-，Zhang X:DEGseq：用于从RNA-seq数据中识别差异表达基因的R包。生物信息学。2010, 26: 136–138. 10.1093/生物信息学/btp612。
第条谷歌学者
Robinson MD，Smyth GK：评估标记丰度差异的中等统计检验。生物信息学。2007, 23 (21): 2881–2887. 10.1093/bioinformatics/btm453。
第条中国科学院谷歌学者
惠特克·L：关于小数泊松定律。生物特征。1914, 10: 36–71. 10.1093/biomet/10.1.36。
第条谷歌学者
Robinson MD、McCarthy DJ、Smyth GK:edgeR：数字基因表达数据差异表达分析的生物导体包。生物信息学。2010, 26: 139–140. 10.1093/bioinformatics/btp616。
第条中国科学院谷歌学者
Robinson MD，Smyth GK：负二项离散度的小样本估计，及其在SAGE数据中的应用。生物统计学。2008, 9: 321–332. 10.1093/生物统计/kxm030。
第条谷歌学者
Cameron AC，Trivedi PK:计数数据的回归分析。1998年，剑桥大学出版社
书谷歌学者
Robinson MD，Oshlack A：RNA-seq数据差异表达分析的标度归一化方法。基因组生物学。2010年11月：R25-10.1186/gb-2010-11-3-R25。
第条谷歌学者
加载器C：局部回归和可能性。1999年，施普林格
谷歌学者
McCullagh P，Nelder JA：广义线性模型。1989年，查普曼和霍尔/CRC，2
书谷歌学者
locfit：局部回归、似然和密度估计。[https://doi.org/cran.r-project.org/web/packages/locfit/]
Agresti A：分类数据分析。2002年，威利，2
书谷歌学者
Engström P，Tommei D，Stricker S，Smith A，Pollard S，Bertone P：使用标记测序对胶质母细胞瘤干细胞系的转录特征进行描述。2010
谷歌学者
Morrissy AS、Morin RD、Delaney A、Zeng T、McDonald H、Jones S、Zhao Y、Hirst M、Marra MA：癌症基因表达谱分析的下一代标签测序。基因组研究2009，19:1825-1835。10.1101/gr.094482.109。
第条中国科学院谷歌学者
Kasowski M、Grubert F、Heffelfinger C、Hariharan M、Asabere A、Waszak SM、Habegger L、Rozowsky J、Shi M、Urban AE、Hong MY、Karczewski KJ、Huber W、Weissman SM、Gerstein MB、Korbel JO、Snyder M：人类之间转录因子结合的变异。科学。2010, 328: 232–235. 10.1126/科学1183621。
第条中国科学院谷歌学者
Benjamini Y，Hochberg Y：控制错误发现率：一种实用且强大的多重测试方法。J Roy Stat Soc B.1995，57:289–300。
谷歌学者
Bullard J，Purdom E，Hansen K，Dudoit S：信使核糖核酸序列实验中标准化和差异表达统计方法的评估。BMC生物信息学。2010, 11: 94-10.1186/1471-2105-11-94.
第条谷歌学者
Bloom JS，Khan Z，Kruglyak L，Singh M，Caudy AA：通过简短测序测量差异基因表达：与双通道基因表达微阵列的定量比较。BMC基因组学。2009, 10: 221-10.1186/1471-2164-10-221.
第条谷歌学者
Smyth GK：Limma：微阵列数据的线性模型。使用R和生物导体的生物信息学和计算生物学解决方案。编辑：绅士R，凯莉V，杜多伊特S，R Irizarry WH。2005年，纽约：施普林格出版社，397-420页。完整文本。
第章谷歌学者
Smyth GK：用于评估微阵列实验中差异表达的线性模型和经验贝叶斯方法。统计应用基因分子生物学。2004年3月：第3条-
第条谷歌学者
Lönnstedt I，速度T：复制微阵列数据。统计正弦。2002, 12: 31–46.
谷歌学者
R：统计计算语言和环境。[https://doi.org/www.R-project.org]
RC、Carey VJ、Bates DM、Bolstad B、Dettling M、Dudoit S、Ellis B、Gautier L、Ge Y、Gentry J、Hornik K、Hothorn T、Huber W、Iacus S、Irizarry R、Leich F、Li C、Maechler M、Rossini AJ、Sawitzki G、Smith C、Smiths G、Tierney L、Yang JYH、Zhang J：生物导体：计算生物学和生物信息学的开放软件开发。基因组生物学。2004年，5:R80-10.1186/gb-2004-5-10-R80。
第条谷歌学者
Bliss CI，Fisher RA：将负二项分布拟合到生物数据。生物计量学。1953, 9: 176–200. 10.2307/3001850.
第条谷歌学者
Clark SJ，Perry JN:通过最大拟似然估计负二项式参数κ。生物计量学。1989年，45:309–316。10.2307/2532055.
第条谷歌学者
Lawless JF：负二项和混合泊松回归。Can J Stat.1987，15:209–225。10.2307/3314912.
第条谷歌学者
Saha K，Paul S：负二项分散参数的偏差修正最大似然估计。生物识别。2005, 61: 179–285. 10.1111/j.0006-341X.2005.030833.x。
第条谷歌学者
快速准确地计算二项式概率。（注意：这是原始论文的副本，不再在线提供。）[https://doi.org/projects.scipy.org/scipy/raw-attachment/ticket/620/loader2000Fast.pdf]
Langmead B、Trapnell C、Pop M、Salzberg SL：短DNA序列与人类基因组的超快和高效记忆比对。基因组生物学。2009年10月：R25-10.1186/gb-2009-10-3-R25。
第条谷歌学者
HTSeq：使用Python分析高通量测序数据。[https://doi.org/www-huber.embl.de/users/anders/HTSeq网站/]
设计。[https://doi.org/www-huber.embl.de/users/anders/DESeq网站]

下载参考资料

致谢

我们感谢Paul Bertone在出版之前共享神经干细胞数据，也感谢Bartek Wilczynski、Ya-Hsin Liu、Nicolas Delhomme和Eileen Furlong共享苍蝇RNA-Seq数据。我们感谢尼古拉斯·德霍姆和朱利安·加涅尔对手稿的有益评论。美国已获得欧盟研究和培训网络“染色质可塑性”的部分资助。

作者信息

作者和附属机构

欧洲分子生物学实验室，Mayerhofstraße 1，69117，德国海德堡
西蒙·安德斯和沃尔夫冈·胡贝尔

作者

西蒙·安德斯
查看作者出版物
您还可以在中搜索此作者公共医学谷歌学者
弗冈·胡贝尔
查看作者出版物
您还可以在中搜索此作者公共医学谷歌学者

通讯作者

与的通信西蒙·安德斯.

其他信息

作者的贡献

SA和WH开发了该方法并编写了手稿。SA实施了该方法并进行了分析。

电子辅助材料

附加文件1：补充。包含所有补充注释和补充图。（PDF 2 MB）

13059_2010_9558_MOESM2_ESM。PDF格式

补充文件2：补编II。PDF文件，以Swave文档的形式显示本文中所有分析的源代码，并附有注释。（PDF 21 MB）

附加文件3：原始数据。Tarball包含所提供分析的原始数据。（TGZ 2 MB）

作者提交的原始图像文件

下面是作者提交的原始图像文件的链接。

图1的作者原始文件

图2的作者原始文件

图3的作者原始文件

图4的作者原始文件

图5的作者原始文件

图6的作者原始文件

图7的作者原始文件

权利和权限

本文根据Creative Commons Attribution 4.0 International License的条款分发(http://creativecommons.org/licenses/by/4.0/)，允许在任何媒体中不受限制地使用、分发和复制，前提是您对原作者和来源给予适当的信任，提供到知识共享许可证的链接，并说明是否进行了更改。知识共享公共领域专用豁免(http://creativecommons.org/publicdomain/zero/1.0/)适用于本文中提供的数据，除非另有说明。

转载和许可

关于本文

引用这篇文章

Anders，S.，Huber，W.序列计数数据的差异表达分析。基因组生物学 11，R106（2010）。https://doi.org/10.1186/gb-2010-11-10-r106

下载引文

收到:2010年4月20日
修订过的:2010年7月22日
认可的:2010年10月27日
出版:2010年10月27日
内政部:https://doi.org/10.1186/gb-2010-11-10-r106