DESeq2的建模原理及简单用法

写在前面的废话

不研究不知道，一研究吓一跳，原来DESeq2这么复杂，这10000多的引用量真不是吹的……

image.png

废话超多系列

DESeq2的差异表达分析涉及多个步骤，具体步骤参见下面流程图中的蓝色部分

image.png

简单地说，DESeq2将对原始reads进行建模，使用标准化因子(scale factor)来解释库深度的差异。然后，DESeq2估计基因的离散度，并缩小这些估计值以生成更准确的离散度估计，从而对reads count进行建模。最后，DESeq2拟合负二项分布的模型，并使用Wald检验或似然比检验进行假设检验。

为什么说DESeq2复杂呢？因为上一篇文章讲了七个步骤，也仅仅只是完成了这个流程图中Estimate size facors这一步。

在使用DESeq2进行基因表达差异分析之前，最重要的是明确我们的研究目的，了解数据中的变异来源。一旦我们了解了数据的主要变异来源，就可以在分析之前提前移除它们，或者通过将这些变量包含在统计模型的公式中对它们进行分析。

DESeq2在进行分析之前需要我们自己书写公式，以便让软件知道变异来源以及在差异表达分析中，我们感兴趣的地方。以下面这个数据为例：

image.png

如果我们知道性别是数据中一个比较显著的变异来源，那么我们就需要将sex写入到统计模型的公式中。公式应该包含数据中的所有因素，这些因素解释了数据中主要的变化来源，其中公式中的最后一个因素，应为我们最为关注的因素

比如，我们想要知道treatment的影响，其中sex和age是主要的变异来源，那么我们的公式则应该为design <- ~sex + age + treatment

公式中的波浪线应该在所有的代数式之前，从而告诉DESeq2在进行差异表达分析时，使用后面的公式。而公式中代数的名称应该与数据框中的列名相匹配。

此外，DESeq2还允许我们研究变异之间的交互作用。比如，我们想知道sex对于treatment的影响，那么我们的公式就应该是design <- ~ sex + age + treatment + sex:treatment

此处需要注意，因为我们关注的是sex对于treatment的交互作用，因此sex:treatment应该放在公式的最末尾

接下来就是无脑运行软件，进行差异表达分析。
首先创建一个DESeq2Dataset对象

dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)

• txi是reads count的矩阵，每一行是一个基因，每一列是一个样本
• colData则是一个因子数据框，每一个因子表示一个样本，相同处理方式的样本采用同样的factor
• design就是刚刚上面所介绍的统计模型的公式

txi和colData的描述可能有点抽象，这里举一个例子进一步说明：
下面就是txi应该有的格式：

image.png

倘若txi如上图所示，则colData则应通过下述代码得到

meta <- factor(rep(c('WT','KO'), each=3))
meta <- data.frame(row.names=colnames(txi), meta)

接下来进行差异表达分析，我们调用DESeq()函数即可

dds <- DESeq(dds)

这一步通过调用DESeq()，将软件的运行结果重新赋给了dds对象。虽然我们仅仅用了一个命令，但是其中涉及到了多个步骤。软件运行的输出信息打印出了它所执行的各个步骤：

estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

除了这种无脑式的一键调用，DESeq2还提供了一些单独的功能，可以让我们一步一步地执行工作流中的每一步。接下来我们详细看看这几个步骤的原理

1. estimating size factors

这一步也就是上一篇文章所说的文库矫正，通过取log，找中位数，减少异常值对scale factor的影响，从而找到一个合适的scale factor。

我们要想单独运行这一步，可以使用函数estimateSizeFactors()，示例如下：

dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)
## 计算sizefactors
dds.sizefactor <- estimateSizeFactors(dds)
## 如果想要知道具体的sizefactor是多少，可以使用sizeFactors()函数
sizeFactors(dds.sizefactor )

sizeFactors()函数除了可以查看表达矩阵评估得到的具体的sizefactor，还可以给一个DESeqDataSet对象的sizefactor赋值，这样DESeq2在对DESeqDataSet对象进行差异表达分析时，就可以使用这个赋值的sizefactor进行后续分析。

dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)
# WT和KO两种处理，六个样本的sizefactor分别是2，2，2，8，9，10
sizefactor <- c(2,2,2,8,9,10)
names(sizefactor) <- c('WT_1','WT_2','WT_3','KO_1','KO_2','KO_3')
sizeFactors(dds) <- sizefactor

所以，如果你的样本加有spike-in，你通过各种方式最终得到了一个scale factor，也可以通过这种方式赋给你的基因表达数据集

2. Estimate gene-wise dispersion

差异表达分析的第二步是对数据离散程度的评估，在RNA-seq的reads count数据中，我们需要知道两点：

1. 为了确定差异表达的基因，我们需要根据组内(重复之间)的方差来确定基因的表达值在组间是否有显著差异
2. 组内(重复之间)的变异需要考虑到方差随表达量的平均值增加的情况，如下图所示(每个黑点是一个基因)。

image.png

为了更加准确的确定差异表达基因，DESeq2需要解释方差和均值的关系。从上图可知，在低表达的基因中，它们的方差也更低，因此DESeq2不希望差异表达基因都是低表达基因。

DESeq2使用离散度(dispersion)作为方差的度量方式，离散度既可以解释基因表达值的方差也可以解释基因的平均表达值。其具体公式为：Var = μ + α*μ^2。其中Var表示方差，μ表示均值，α表示离散度。因此我们可以得到这么一个关系

	离散度
方差增加	离散度增加
平均值增加	离散度降低

那么在表达水平较高的基因中，离散度的平方根 就等于方差系数。其中σ是标准差，μ是平均值。那么α=0.01就意味着，在样本生物重复之间存在着10%的标准差。

表达水平较高时，μ对于公式的影响显著小于

因此，具有相同平均值的基因的离散度仅根据其方差而存在差异，离散度反映了一个给定平均值的基因表达的差异程度。

那么接下来一个比较重要的问题，就是如何将离散度与我们的模型建立联系呢？为了更精确的为我们的数据建模，我们需要知道每个基因组内方差的精确评估值。

但是，生物这个行当，一般3个重复就了不起了，6个重复就顶破天了。看起来似乎挺多的，但是远远不够，这就导致我们得到的组内方差是相当的不可靠……

6个重复都不一定够，所以那些做RNA-seq一个处理只有一个重复的同学，你们是想搞事情么？

image.png

针对这个问题，DESeq2使用一种叫作shrinkage的方法，共享基因之间的表达信息，根据基因的表达水平生成更为准确的方差估计。DESeq2假设具有相似表达水平的基因也具有相似的方差。这样DESeq2就可以基于基因的平均表达水平和离散度来评估每个基因的离散度。

在下文里我就暂且把shrinkage翻译为收缩

3. Fit curve to gene-wise dispersion estimates

DESeq2的第三步就是根据基因的离散度拟合一个曲线。那么为什么要做拟合呢？不同的基因生物学重复中存在不同的方差，但是，在所有的基因中，将会有一个合理的离散分布。

这个曲线如下图红线所示，其中红线的横坐标是基因的表达强度，纵坐标是理论离散值。而每一个黑点的横坐标是基因的平均表达水平，纵坐标是经过最大似然估计的离散值。

image.png

4. Shrink gene-wise dispersion estimates toward the values predicted by the curve

有了拟合曲线，接下来就是对基因表达水平的离散度进行矫正，即：将基因的实际离散度向红线收缩(shrink)。当样本量较小时，该曲线可以让我们更为准确的识别差异表达基因。既然知道要将基因的离散度向红色曲线收缩，那么收缩多少比较合适呢？有两点需要考虑：

1. 基因的离散度距离红色曲线的距离
2. 样本量（样本量越大，则收缩的越少）

这种方法在差异表达分析时，可以极大的减少数据的假阳性。离散度较低的基因朝着理论值收缩，从而得到一个更为准确的离散值。而那些离散度较高的基因，则不能无脑朝着理论值收缩。

这是因为这类基因可能不遵循建模假设，并且由于生物学或技术原因使得这类基因与其他基因具有更高的方差。DESeq2识别这类基因后，将不采用shrink方法对它们进行处理。下图中蓝色圆圈圈住的便是这类基因。

image.png

5. GLM fit for each gene

DESeq2的第五步，对每个基因使用广义线性模型进行拟合，我自己也没搞明白，就不在这里献丑了

写在后面的话

也许你会问既然DESeq()函数可以一键实现这么多操作，为什么还要一步一步单拎出来做呢？主要是我自己有点强迫症，知其然不知其所以然的感觉太痛苦了……更何况，有了spike-in的教训，我还敢无脑运行软件么？

参考资料

• https://hbctraining.github.io/ ：Introduction to DGE
• Beginner's guide to using the DESeq2 package