基因芯片数据分析（六）：DESeq2包的基本原理

文库均一化面临的问题

DESeq2是另外一个分析差异基因的R包，它的功能很多，使用也比较复杂。我们在前面提到过，RPKM，FPKM与TPM是常用的用于均一化不同的样本reads数的方法，不过DESeq2和edgeR并不使用前面的三种方法，因为在对文库进行均一化时，存在两个问题，如下所示：

第1个问题就是，不同样本的文库大小还有差异，需要对其进行调整。我们先看一个简单的案例，在这个案例中，我们假设某个物种的基因组上只有6个基因，如下所示：

现在我们手中有这个物种的两个测序样本，其中Sample #1的reads数为635，Sample #2的reads数为1270，它们的各自基因对应的reads数如下所示，从中我们可以看到，Sample #1的reads总数是Sample #2的reads总数的一半，如下所示：

我们再看一下这两个样本每个基因对应的reads数，从中我们可以发现，Sample #2中每个基因对应的reads数也基本上是Sample #1中每个基因对应的reads数的2倍，这种差异不可能是生物学的原因，因为物种都一样，出现这种现象的原因就是测序深度，而RPKM，FPKM和TPM都能处理这样的问题（也就是测序深度不同造成的reads不同），如下所示：

第2个问题是，我们需要调整由于文库补偿（library composition没有找到相应的中文翻译，此处译为“文库补偿”）造成的差异。RNA-seq或者是其它的高通量测序技术通常会比较不同组织类型之间的测序数据差异，例如我们可能会比较肝脏与脾脏的差异。这个时候就可能会出现问题，因为在肝脏中，存在着某些特异性基因，这些基因只在肝脏中大量转录，而在脾脏中转录不活跃，这就是文库补偿造成的差异的一个案例，如下所示：

我们看一下A2M这个基因，在Sample #1中，这个基因的reads数是635，而在Sample #2中这个基因的reads数是0，这就说明，A2M这个基因只在Sample #1中转录，它的reads数是563，但是在Sample #2中，这个563的reads数必定是分配到其它的基因上去了，如下所示：

现在我们看一下Sample #2中所有基因的reads数，我们可以发现，在Sample #2中，A2M的reads数为0，如下所示：

如果出现了这种情况，那么RPKM，FPKM和TPM是无法处理的，此时就需要使用DESeq2（或者是edgeR）来处理，DESeq2能处理各种测序的数据集，它能解决的两个问题：

第一，不同文库大小之间的差异；

第二，文库补偿效应造成的差异。

如下所示：

DESeq2均一化的步骤

现在我们看一个简单的案例，在下面的这个案例中，我们有3个样本，每个样本有3个基因，在这个案例中，我们的目标是计算每个样本的标准化因子（scaling factor），标准化因子会解决测序深度（read depth）和文库补偿（library composition）的问题，如下所示：

第一步：对reads数取自然对数

DESeq2中默认是对reads数取以e为底的对数（但是也可能设置以2或10为底的对数）如下所示：

自然对数比较好理解，以Sample #2的Gene 1为例，它的reads数是10，那么loge(10)=2.3，另外需要注意的是，如查reads数为0，那么对数就是-Infinity，缩写为-Inf，表示负无穷，如下所示：

第二步：求所有样本中，相同基因对数的均值

在这一步中，我们需要对同一个基因在所有样本中的数值取均数，以Gene 1为例，Gene 1在这三个样本中的数值分别为-Inf，2.3，1.4，由于-Inf是负无穷，因此加起来也是-Inf，再除以3，还是这个数，我们再看Gene 2的数值，它在3个样本中的数值分别为0.7，1.8，2.5，那么它的均值就是1.7，如下所示：

用对数的一个原因就是对数不太容易受到异常值的影响，我们看一下Gene3，它在Sample #3中的原始reads数是200，而在Sample #1和Sample #2中的数值分别为33和55，因此200明显是一个异常值，如果直接使用reads数来计算均数，那么结果就是96，我们再看一下使用对数的均值，计算结果是4.3，这个4.3是指数，现在把它转换成与reads数对应的原始数值，也就是e的4.3次方，结果为73.7，与96相比，前者还是比较小的，说明使用对数的受异常值的影响比较小（这种先取对数，后用对数求均值的方法叫几何平均数），如下所示：

第三步：去除掉Infinity

在这一步中，我们要把在第二步中计算出的含有Inf的基因给剔除掉，也就是说要在样本中把reads数为0的基因给剔除掉，如果如果我们比较的是肝脏和脾脏的转录组，按这种剔除方法，我们会把那些所有只在肝脏（或脾脏）中转录的基因都给剔除掉，从理论上讲，最终剩下的基因基本上就是管家基因（house keeping genes）了，也就是说在不同组织类型中都表达相似的基因，如下所示：

第四步：矩阵减均值

在上一步中，我们把数据中的Inf值去除了，在这一步中，把每个样本中的每个基本减去该基因在所有样本中的均值，以Gene 2为例，它在经对数转换后的Sample #1、Sample #2、Sample #3中的数值分别为0.7，1.8和2.5，它在所有样本中的均值为1.7，这3个数值都减去1.7，就是-1.0，0.1，0.5，其余基因操作均如此，如下所示：

在这一步骤中，我们用的是对数转换后的数值相减，其本质上是以某个基因（这里称为基因X）的平均值为参考，对每个样本中的基因X进行均一化，就是原始reads数的相除，如下所示：

第五步：计算每个样本的中位数

在这一步中，我们要计算每个样本中，所有基因的对数的中位数，以Sample #1为例，在这个样本中，它所有的基因只有Gene2和Gene3，分别为-1.0和-0.8，它的中位数就是-0.9，其余的样本的计算也是如此，这里需要注明一下，因为在这个案例中，一个样本的所有基因只有2个，因此它的中位数和无数是相等的，但是如果基因的数目很大，这两个值就不一定相等了，例如某批数据，分别为1，3，5，7，10000，那么这批数据的中位数就是5，均值则为2003.2，这个均值受异常值的影响很大，再看一批数据，分别为-1000，3，5，7，9，那么这批数据的中位数还是5，均值就变成了-195.2了，从中我们可以发现，均值比较容易受异常值的影响，但中位数对异常值则不敏感。

这里使用中位数主要是为了排除一些极端表达基因的影响，极端表达指的是表达量极高或极低，它们能够对均值造成很大的影响，而那些表达量差异极大的基因对于中位数的影响，并不比那些表达量差异较小的基因对中位数的影响大。再加上那些表达量差异极大的基因数量通常情况下很少，因此，更多情况下，中等程度表达差异的基因和管家基因对中位数的影响更大一些，如下所示：

第六步：将中位数转换为真数，计算每个样本最终的标准化因子

在这一步中，我们要把第五步中的对数中位数转换为真数，如下所示：

此时，我们就得到了这3个样本的标准化因子了，此时进入第七步。

第七步：原始reads数除以标准化因子

在这一步中，我们要把原始的reads数除以这个标准化因子，以Sample #1为例，它的标准化因子为0.4，那么Gene 1，Gene2和Gene 3的reads数分别为0，2，33，它们分别除以0.4，结果为0，5，82.5，由于read不可能是小数，四舍五入，结果就是0，5，83，其余的样本均按此方法处理，如下所示：

DESeq2文库标准化因子总结

log转换会剔除那些只在某个样本类型中的转录的所有基因（例如肝脏与脾脏），这种处理也会消除原始reads数的异常值（通过几何均数）。

中位数的处理会进一步降低那些高数值reads数基因的影响，从而关注那些中等表达程度的基因，如下所示：

不相关过滤（Independent Filtering）

在这一部分中，我们会讲一下对那些低reads数的基因过滤的方法，这个过程又叫不相关过滤（Independent Filtering没有找到相应的中文译名，这里译为“不相关过滤”），也叫解决多重检测问题，如下所示：

有关假阳性

当我们做一个统计学检验时，就有一定的概率得到的是错误结论，简单来说就是，我们认为“p值小于0.05就是有统计学差异的”，也可以这么说，“我们有5%的概率会得到一个假阳性的结果”，如下所示：

我们来用一个简单的案例来说明这个问题，在下面的这个案例中，其中红点表示突变组，黑点表示对照组，我们检测了Gene 1和Gene 2的这两个基因，它们对应的p值如下所示：

从上面的结果可以看出，它们的p值分别为0.03和0.13，基因数目是2，2的5%是一个很小的数值，因此，我们很不太可能得到一个假阳性的结果，但是，如果我们比较了整个基因组中的所有基因（大概有2万个），用于观察哪些基因调控出现了异常时，情况就不太一样了，如下所示：

2万的5%就是1000，也就是说，如果我们以p值0.05为阈值进行统计的话，我们有可能得到约1000个假阳性结果，对于这个问题，我们通常会采用FDR和Benjamini-Hochberg方法来解决，不过，即使是这样，我们还会遇到另外的一个问题，先看一个案例。

在下面的图形中，我们看到了两个独立的分布，其中红色曲线表示X品系小鼠的体重，蓝色的曲线表示Y品系小鼠的体重分布，如下所示：

如果我们称量一个X品系小鼠的体重，那么它有很大的概率会落到红色曲线的中间部分，如下所示：

如果我们称量3只X品系的小鼠，它们都有可能落在红色曲线的中间部分，如下所示：

同样类似的，还有Y品系小鼠的体重，如果我们称量3只Y品系的小鼠，那么它们也有很大的概率会落在蓝色曲线的中间区域，如下所示：

如果我们对这两批数据进行t检验，得到的p值小于0.05的话，我们就可以正确地下一个结论，即这两批数据存在着差异，这种差异是来源于两种小鼠的体重分布不同，也就是说这是由于小鼠品系不同千万的体重差异，如下所示：

但是，有时候我们称量的小鼠体重会出现重叠，此时，我们进行t检验，它的p值有可能就是大于0.05的，这就是假阴性，如下所示：

如果我们使用计算机从这些分布中进行1000次抽样（一次抽样包括从对照组中抽3只，突变组中抽3只），如下所示：

现在做1000次t检验，我们看一下这1000次的p值，如下所示：

其中我们可以发现，有949个结果的p值是小于0.05的，有51个假阴性（p值大于0.05），如下所示：

目前为止，我们进行的每次检验都被认为应该是一个真阳性，它的p值是小于0.05的，这句话的意思是说，我们希望检验的这两个数据（也就是突变组与对照组的数据）是有差异，当然了，实际结果可能没有差异，这就是假阴性。出现真阳性是因为我们每次检验使用的样本都是来源于不同的分布，如下所示：

当我们加入一些没有统计学差异的样本时，也就是说来源于同一个分布的样本时，它的p值就会大于0.05，如下所示：

但是，即使是来源于同一分布的两组数据，有的时候，也会出现p值小于0.05的情况，就像下面的这个样本，如下所示：

现在我们绘制一个直方图，这个直方图是2000个p值，其中1000个来源于两个不同的分布的两组数据的t检验的p值，另外的1000个是来源于同一个分布的两组数据的t检验的p值，如下所示：

从上面的直方图我们可以知道，其中有993个小于0.05的p值，如下所示：

在这993个小于0.05的p值中，有949个是真阳性（这是由第一组p值构成的，也就是来源于两个不同分布的两组数据的t检验的p值），有44个假阳性的p值（这是由第二组p值构成的，也就是来源于同一个分布的两组数据的t检验的p值），如下所示：

由于在这些小于0.05的p值中，只有约4%的p值是假阳性（44/993=4.4%），因此我们并不需要FDR校正，但是，如果在实际的分析中，我们并不知道p值的构成，因此就需要进行FDR校正，如下所示：

我们来看一下，这个数据经过FDR校正后的结果，校正后，有846个p值仍然小于0.05，在这846个数值中，有827个是来源于原来真阳性的949个，所占比例为89%，有19个是假阳性，占846的比例为2%，如下所示：

现在我们把数据的构成改变一下，让它更像RNA-seq的数据，此时我们增加样本的数目到6000个，其中1000个来源于两个不同的分布，5000个来源于相同的分布，这5000个样本进行t检验，它们的p值有更大的可能性大于0.05，如下所示：

绘图出样本t检验的p值直方图，如下所示：

现在我们看一下这些p值的分布，其中有1215个p值小于0.05，其中949个p值是真阳性，有266个p值是假阳性，假阳性的比例为22%，此时就需要进行FDR校正，如下所示：

经过FDR校正，我们发现，只有256个校正后的p值小于0.05，其中250个来源于原来的949个真阳性，其比例为26%，剩下的6个是假阳性，来源于原来256个假阳性，比例为2%，从这个结果我们可以看出来，FDR能够限制假阳性的数目，但是它同时还能减少真阳性的数目，如下所示：

现在我们把样本数目增加到11000，其中1000个来源于两个不同的分布的两组数据，10000个来源于相同的分布的两组数据，现在对这些样本进行t检验，将它们的p值绘制成直方图，如下所示：

其中1430个p值小于0.05，跟前面的结果一样，有949个真阳性的p值，但是有481个假阳性的p值，假阳性的比例高达34%（481/1430），如下所示：

经过FDR校正，只有56个p值仍然是小于0.05，其中54个真阳性的p值，2个假阳性的p值，也就是说，原来949个真阳性中，只剩下了54个真阳性，只有原来的6%，而假阳性占到FDR校正后的小于0.05的p值数目的4%，从前面的计算结果可以看出来（2000个样本，6000个样本，11000个样本），每当我们增加样本的检验数目时，通过FDR校正的真阳性的（p小于0.05）数目都会减少，如下所示：

用折线图表示就是下面的这个样子，其中绿色曲线表示的是通过FDR校正后，p值仍然小于0.05的数目，橘黄色的线表示经过FDR校正后，p值仍然小于0.05的假阳性的数目，这张图表说明，无论检验的数目是多少，FDR总能控制假阳性的比例，这在我们的模拟数据中已经得到了很好验证，不过随着检验数目的增多，真阳性的比例却在下降，如下所示：

这也说明了，Benjamni-Hochberg方法时行的FDR校正还有很大的改进空间，edgeR和DESeq2都有自己的方法来对前面的模拟检验进行过滤。

edgeR和DESeq2的过滤

edgeR和DESeq2对检验结果进行过滤的基本思路就是，那些reads数极少的基因提供的信息量有限，可以把它们从数据集中剔除掉。换句话讲，如果这个基因在一个样本中只有1个或2个转录本，在另外一个同样的样本中只有3或4个转录本，那么就很难得到它们的精确reads数，即使这些低reads数的基因有着生物学意义，那么我们也能剔除它们，如下所示：

edgeR的过滤

在了解edgeR过滤数据之前，我们先了解另外一个概念，即CPM，CPM的全称是Counts Per Million，它用于校正不同文库之间的由于测序深度造成的差异，edgeR是按照CPM来剔除数据的，它的阈值是1，edgeR会保留那些在2个及2个以上的样本中CPM大于1的基因，剩余的基因就被剔除掉，具体的后面我们会讲到，如下所示：

CPM

下图说明了如何计算CPM，流程为：①总reads数除以100万；②用每个样本中的每个基因对应的reads数再除以①中的数字，就是CPM。

总reads数除以100万主要是为了方便计算CPM，否则CPM的数值会非常小，不方便，如下所示：

现在我们有了所有样本中的所有基因的CPM值，此时就要保留那些有2个或2个以上样本中，CPM值大于1的基因，剩下的基因除去，如下所示：

现在我们看一下如何保留需要的基因，一个基因一个基因地来看：

gene 1，它在kidney 1和kidney 2这两个样本中的CPM值都大于1，因此保留，如下所示：

再看gene 2，gene2在kidney 1和kidney 2这两个样本中，CPM值都大于1（有一个是等于1），因此我们保留基因2，如下所示：

再看gene 3，由于gene 3只在一个样本，即liver 2中的CPM值大于1，即使这个值很大，我们也不需要（下图中的这个值是103.3，跟前面视频的数据不一样的，我估计作者在这里主要是为了说明一个问题，即使某个值在一个样本中很大，在其他样本中仍然小于1，也是要剔除这个数据），如下所示：

再看Gene 4，由于有2个样本的CPM值小于1（可能作者没有算另外两个CPM等于0的样本），因此也要剔除，如下所示：

再看Gene 5，它在所有的4个样本中CPM值都小于1，因此Gene 5也不用保留，如下所示：

再看Gene 6，这个基因在2个样本中CPM值大于1，即使这两个样本不是同一类型的，但我们可以保留这个基因，如下所示：

edgeR的这种方法很简单，但我们要知道，测序深度会影响这个数值，例如如果一个样本中有500万条reads，那么CPm的标准化因子就是500万/100万=5，如果我们有5条reads比对到某个基因上了，那么我们就会知道这个基因是1CPM，即5/5=1CPM，如果我们有800万个基因，那么标准化因子就是80，1CPM=80 reads，这个数值就非常大了，如下所示：

换句话讲，有时候，我们也需要一个比较大的CPM阈值了，例如，如果我们的一个样本中有5万条reads（单细胞测序的reads数可能就这么多了），那么它的CPM标准化因子为0.05，即5万/100万=0.05，如果一个基因有一条reads比对上了，那么这个基因的CPM就是20，即1/0.05=20。即使一个基因在生物学水平进行了转录，但是你只能低到一个低的reads数和低的总reads数，它仍然存在着很大的噪音，如下所示：

edgeR中CPM阈值的确认

此时，我们可能面临一个问题：如何找一个很好的阈值（cutoff）？

我们来看一个真实的案例。

这个案例是从别人那里得到的一批真实数据，在这批数据中，每个样本有2200万条reads（8个样本，其中包括4个野生型，4个敲除型）；

在不剔除任何基因的情况下，运行edgeR，得到一批原始的p值；

使用不同的CPM阈值（cutoff）来过滤基因，然后对这些p值进行校正，如下所示：

此时我们使用不同的CPM阈值来过滤一下基因，然后再校正p值，就是下图所示的内容，其中y轴表示经过FDR校正后，p值仍然小于0.05的基因的数目，x轴表示最小的CPM的阈值，如下所示：

其中横坐标上的0表示，没有剔除掉任何一个基因，如下所示：

如果我们选择1，也就是Minimum CPM Threshold的值为1（也就是推荐的值），此时FDR小于0.05的基因的数目大概是3400个，在此处的reads数还是比较多的，这说明这个阈值还是比较严格的，如下所示：

如果我们使用一个较低的阈值，我们就会发现，此时FDR小于0.05的值基因数目大概是3600个，比CPM=1时的数目多了约200个，如下所示：

总之，edgeR的主旨就是：仔细选择阈值（cutoff），你需要在计算了p值之后再尝试不同的CPM阈值，如下所示：

DESeq2中CPM阈值的确认

在确认CPM的阈值方面，edgeR和DESeq2有两点不同。

第一，edgeR计算CPM的阈值时，它会检查每个基因在不同样本中的数值，从而确定保留下来的基因中，至少有2个样本的CPM值大于阈值，在下面的这个案例中，我们就发现，gene 1在kidney 1和kidney 2这2个样本中，它的CPM值都大于1，因此保留gene 1，如下所示：

相比之下，DESeq2则是要检查某个基因在所有样本中均一化后的数值的均值，还是以上面的案例说明，如果是DESeq2，它比较的结果就是，这4个样本的平均值都在阈值之上，因此可以保留gene 1，如下所示：

此时，你或许会想，如果我使用了DESeq2的方法，那么对于异常值，该如何处理？

看下面的案例，gene 1在kidney 1，kidney 2，liver 1中的CPM都是0，而在liver 2中则是5000，平均值是1250 CPM，按照edgeR的标准，gene 1不会被保留，因此它只有一个样本的CPM大于1，如下所示：

但是对于DESeq2来说，gene 1在4个样本中的平均值为1250，它会被DESeq2保留（这种情况下DESeq2无法处理，也就是一个类型中有2个样本），如下所示：

但DESeq2本身具备异常值的检测方法，不过这种方法只有在每类别中，有超过2个样本的时候才发挥作用（具体的可以参考这篇文献：Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2）。

下图是分别使用这两种方法（edgeR和DESeq2）对同一批数据进行计算的结果，其中横坐标表示最小CPM阈值，纵坐标表示有显著差异的基因的数目，如下所示：

从中我们可以发现，这两种方法基本上都有在同一区域达到顶点，如下所示：

这两种方法同时会产生一个相似的阈值，如下所示：

现在我们看一下这两种方法不同之处，其中一处区别在于x轴，如下所示：

再看另外一处差异，DESeq2能绘制X轴是百分位数（用于替代min CPM 阈值），Y轴是显著差异基因数目的二维坐标系，如下所示：

根据这种二维坐标图，我们可以看出来，0表示在这个阈值之下的基因是0%，如下所示：

0.2表示在这个阈值之下，有20%的基因，如下所示：

分位数（Quantiles）非常有用，因为我们可以看到，CPMs的计算依赖于测序深度。无论这个文库是800万条reads，还是8000万条reads，10%的基因总是总是小于0.1百分倍数，如下所示：

我们可以选择使用分位数法和最小CPM，如下所示：

现在看第三处差异，DeSeq2会对这些点进行拟合，找到一曲线，如下所示：

DESeq2就会找到这条拟合曲线的最大位置处，如下所示：

阈值就在这条曲线的最大位置处，关送去拟合曲线和原始值之间的标准差，换句话讲，在峰值的噪声范围内，第一个分位数就是CPM的阈值，如下所示：

如果在这个阈值之上没有原始值，那么就不过滤，如下所示：

至此，我们就知道了edgeR和DESeq2是如何过滤基因的，其中edgeR是保留那些在2个或2个以上样本中，CPM大于最小阈值的基因，而DESeq2则是，保留那些平均CPM大于最小CPM的基因，然后绘制显著基因与分位数的散点图，找到拟合曲线，再用最大值减去噪声，即是阈值，如下所示：

不过现在我们知道，我们可以同时联合使用这两种方法，如下所示：

建议

如果你使用edgeR，㙘在计算出了p值之后，找到CPM的阈值。在edgeR的基因选择标准方面，使用DESeq2可以很容易地找到找到最优的CPM阈值。

如果你使用DESeq2时，如果在每个分类（例如野生型组和突变型组）只有2个样本，此时要尤其注意异常值，如下所示：

好了，关于DESeq2和edgeR包，我们上一篇文章和本文已经介绍，后面我们将介绍limma包，以及这些包的实战教程，扫码关注我们，一起学习！同时也欢迎批评指正！！！