Differential expression analysis for sequence count data

摘要

高通量测序分析以计数数据的形式提供定量读数。为了正确推断此类数据中的差分信号并具有良好的统计能力，需要估计整个动态范围内的数据可变性和合适的误差模型。文章提出了一种基于负二项分布的、方差和均值通过局部回归联系起来的方法。

研究背景

DNA 片段的高通量测序用于一系列定量分析。这些检测的一个共同特征是它们对大量 DNA 片段进行测序，这些片段反映了例如生物系统的 RNA 分子库（RNA-Seq [ 1 , 2 ]）或核苷酸结合分子的 DNA 或 RNA 相互作用区域。通常，这些reads会根据它们与目标基因组的共同区域的映射分配给一个类别，其中，对于RNA-Seq，每个类代表一个目标转录物，对于ChIP-Seq，每个类代表一个结合区域。
我们希望使用统计测试来确定，对于给定的基因，观察到的reads计数差异是否显著，也就是说，它是否大于仅仅由于自然随机变异而产生的预期值。

reads

如果从具有给定固定基因部分的群体中独立抽样reads，则reads计数将遵循多项式分布，该分布可近似为泊松分布。
现在泊松分布已经用于基于测序表达中。泊松分布只有一个参数，由其平均值唯一确定；它的变化和所有其他属性都来自它；特别是，方差等于平均值。
但在实际使用中，我们发现泊松分布的假设过于严格：它预测的变化比数据中看到的要小。由此产生的统计测试不能控制所宣传的I型错误（错误发现的概率）。
为了解决这个所谓的过度分散问题，有人提出用负二项式 (NB) 分布对计数数据进行建模。
NB 分布具有由均值μ和方差σ 2唯一确定的参数。
感兴趣的数据集中的重复次数往往太少，无法可靠地估计每个基因的平均值和方差这两个参数。
因此我们可以假设平均值和方差通过s2=μ+aμ2进行关联，其中一个比例常数a在整个实验过程中是相同的，并且可以从数据中进行估计。因此，每个基因只需估计一个参数，就可以把这个分布应用于少量重复的实验。

reads数量建模

假设分配给基因i 的样本j中的reads数量可以通过负二项式 (NB) 分布建模：

它有两个参数，均值和方差。读取计数$K_{ij}$是非负整数。
在实验中，我们不知道${\mu }_{ij}$和${\sigma }_{ij}^2$。我们需要从数据中估计它们。
但有用的数据重复次数很少，需要进行进一步的建模假设以获得有用的估计。
因此文章提出了基于一下三个假设的方法。
一、
首先，平均数${\mu }_{ij}$，即样本j中基因i的观察计数的期望值，是依赖于条件的每个基因值$q_i$，${\rho }_{(j)}$（其中${\rho }_{(j)}$是样本j的实验条件）和尺寸因子$s_j$的乘积：

$q_i$,${\rho }_{(j)}$与条件${\rho }_{(j)}$下基因i片段的真实（但未知）浓度的期望值成正比。size factor$s_j$表示样本j的覆盖范围或采样深度，我们将使用术语common scale表示数量，例如q i , ρ ( j )，通过除以s j调整覆盖范围。
二、
方差${\sigma }_{ij}^2$是散粒噪声项和原始方差项的总和。

三、
假设每个基因的原始方差参数$v_{i,\rho }$是$q_{i,\rho }$的平滑函数。
这个假设是必要的，因为重复的数量通常太少，无法仅从该基因的可用数据中精确估计基因i的方差。这个假设允许我们将来自具有相似表达强度的基因的数据汇集起来以进行方差估计。

模型拟合数据

数据是一个n × m的计数表$k_{ij}$，其中i = 1,…, n，对基因进行索引，而j = 1,…, m对样本进行索引。该模型具有三组参数：

1. m 个尺寸因子$s_j$，来自样本j的所有计数的期望值与$s_j$成正比。
2. 对于每个实验条件$\rho ,n$表达强度参数$q_{i\rho }$；它们反映了条件$\rho$下基因i片段的预期丰度，即基因i计数的期望值与$q_{i\rho }$成正比。
3. 平滑函数$v_{\rho }:R^{+}\to R^{+}$; 对于每个条件$\rho ,v_{\rho }$对原始方差$v_{i\rho }$对预期均值$q_{i\rho }$的依赖性进行建模。

差异表达测试

假设有$m_A$用于生物条件A 的重复样本，和$m_B$个用于条件B 的样本。
对于每个基因i，想权衡该基因在两种条件之间差异表达在数据中的证据。
也就是说要检验空假设：$q_{iA}=q_{iB}$。
其中$q_{iA}$是条件 A 样本的表达强度参数，$q_{iB}$是条件 B的表达强度参数。为此，我们将每种情况下的总计数定义为检验统计量:

他们的总和为：$K_{iS}=K_{iA}+K_{iB}$。
在零假设下，我们可以计算任意一对数字a和b的$K_{iA}=a$和$K_{iB}=b$事件的概率。可以用$p(a,b)$表示这个概率。
一对观测计数和$(k_{iA},k_{iB})$的P值为小于或等于$P（k_{iA}，k_{iB}）$的所有概率之和，假设总和为$k_{iS}$：

上述总和中的变量a和b取值 为$0,...,k_{iS}$。

$p(a,b)$如何计算呢？
首先，假设在零假设下，来自不同样本的计数是独立的。然后，p $p(a,b)=Pr(K_{iA}=a)Pr(K_{iB}=b)$。
因此，问题就转换成计算事件$K_{iA}=a$的概率，和$K_{iB}=b$的概率。随机变量$K_{iA}$是$m_A$ 负二项分布随机变量的总和。

我们通过NB分布近似计算其分布，其参数从$K_{ij}的参数中获得。为此，首先根据两种情况的计数计算合并平均估计:

这样就解释了空假设规定$q_{iA}=q_{iB}$的事实。条件A的总平均值和方差为：

使用的数据集

果蝇胚胎中的 RNA-Seq：在这个数据集的每个样本中，一个基因被设计为过度表达，我们比较了两个这样的条件中的每一个的两个生物学 replicates。表示为“A”和“B”。
神经干细胞的标签序列。
酵母的 RNA-Seq。
HapMap 样本的 ChIP-Seq。

方差估计

主要是验证了DESeq保持对 I 类错误的控制。
我们将果蝇数据中条件A的一个 replicates与另一个进行对比，对两个样本使用从两个 replicates估计的方差函数，绘制了变异系数 (SCV)，即方差与均方的比值。

图1

图1中：
橙色线是拟合的$w（q）$。
紫色线显示了泊松分布对两个样本中每一个样本的隐含方差。
橙色虚线是edgeR使用的方差估计。
(b) 与 (a) 中的数据相同，重新调整y轴以显示变异系数 (SCV) 的平方，即所有数量均除以均值的平方。