StatQuest生物统计学专题 - library normalization进阶之DESeq2的标准化方法

RPKM,FPKM,TPM等标准化方法还有那些问题？DESeq2的标准化方法的原理就是提高中等表达基因的地位一个例子

在上一节的StatQuest生物统计学专题中，我们简单直白的讨论了RPKM,FPKM,TPM的定义和生物学意义，明白了RPKM,FPKM,TPM标准化方法就是为了去除基因长度和测序深度对测序Read数的影响，见StatQuest生物统计学专题 - RPKM,FPKM,TPM，并且在一般意义上，更推荐使用TPM标准化方法。

然而不得不说明的是，诸如DESeq2、edgeR等差异表达分析软件都不是使用的RPKM,FPKM,TPM方法，为什么呢？

RPKM,FPKM,TPM等标准化方法还有那些问题？

其实TPM之类的标准化方法虽然解决了基因长度和测序深度的影响的问题，但还是不能解决一个问题：那就是测序文库组成不同造成的差异。

什么意思呢？

我们知道RPKM,FPKM,TPM是可以解决由于测序深度的差异而引起的Read数变异，如下图所示，样本2#的总Read数是样本1#的2倍，每个基因的Read数也是1#的2倍。我们知道这种Read数差异并不是基因表达的不同，而是由于测序深度不同所致，只需要将样本1#、2#除以各自的总Read数，那么这个Read数变异就会得到修正。

让我们再考虑两种新的情况，我们知道RNA-seq或其他高通量方法，往往是不同组织间的基因表达比较，比如肝细胞或脾细胞，他们会有一些各自特异性表达的基因，或者对于一些敲减基因的实验，有些基因在其中一个样本中高表达，而在另一个样本中不表达。

如下图所示的两个样本，两者的测序深度是一样的，总Read数相同。但是由于A2M基因在样本2#中不表达，导致563个Read会分配到样本2#的其他5个基因中去：如样本2#的基因A1BG的Read数是235，大于样本1#中的30。然而实际上这种差异并不是生物学效应所致，而是由于基因A2M被敲减所致，并不是A1BG等5个基因在样本2#中表达增加了，而是基因A2M在样本2#中表达减少了。

在这种情况下，这种测序文库组成不同的差异是RPKM,FPKM,TPM等方法无法解决的。

DESeq2的标准化方法的原理就是提高中等表达基因的地位

那么如何解决这个问题呢？

本周先看一下DESeq2是如何进行library normalization的，DESeq2的标准化方法共有7步，看起来很繁琐，但是原理很简单，它有一个贯穿始终的基本思想——提高中等表达基因的地位。

而且这7步只是为了得到一个标准化因子，并进行变换。

首先以下述数据集为例，共有3个样本，每个样本有3个基因：

第一步 对Read矩阵取对数变换

DESeq2默认是使用自然对数，也可以使用log2或log10。

第二歩 取各基因的平均数

要说明的是，由于样本1的基因1的Read数是0，所以在取对数时，它的值是-Inf（负无穷大），因此对基因1取平均数时就直接得出是-Inf即可。

为何要取对数？

其实是为了减少高表达异常值对标准化的影响，需要注意的是异常值不代表是错误值，只是说它相比数据趋势比较异常。

以原始表达矩阵为例，基因3的3个Read数分别是33、55、200，尤其是200，相比较整个表达矩阵来说是一个高表达异常值。

如果对原始表达矩阵求基因均值，那么基因3的均值是(33+55+200)/3 = 96，而对于对数变换后的表达矩阵来说，基因3的均值是(3.5+4.0+5.3)/3 = 4.3； e^4.3 =73.7，73.7<96，也就是说对数变换后基因3的均数受到200的影响更小（这种取对数求得的平均数是几何平均数）。

第三步 过滤掉-Inf基因

将存在-Inf值的基因过滤掉，过滤掉的基因不再参与标准化因子的计算。

实际上，这一步是把在一个或多个样本中存在零表达的基因剔除。假定本实验是在比较不同组织细胞如肝细胞和脾细胞的表达量差异，那么这一步会剔除掉组织特异性表达的基因，而只保留管家基因——在不同细胞中都或多或少会表达的基因。

第四步 将对数矩阵减去对数均值，得到对数比值矩阵

将对数矩阵的每个Read分别减去此基因的对数矩阵。

对数相减的意义何在？

对数相减其实是真数相除，也就是：

log(reads for gene X)-log(average for gene X) = log(reads for gene X/average for gene X)

注意：此时的average for gene X是几何平均数。

第五步 计算每个样本的对数比值矩阵的中位数

取中位数，而不是均值，也是为了进一步降低异常值的影响，具有较大表达差异的基因对中位数的影响甚微。考虑到绝大部分情况下，表达差异大的基因都是很少的，所以这个“中位数”更能代表的是中等表达基因或管家基因的情况。

第六步 将对数中位数转换为其相应的真数，得到各个样本的标准化因子

第七步 将原始表达矩阵除以这个标准化因子

原始矩阵的每个样本的全部Read数均除以各自的标准化因子（包括较大表达差异的基因）。

我们可以看到标准化之后，对于基因3来说，样本1#的Read数提升，而样本3#的Read数下降，基因3在3个样本中的表达其实是接近的。

做一下总结：

1. 对数变换可以减少只表达在某些样本中的基因的影响，同时还可以减少异常值的影响（几何平均数）；

2. 将在某些样本中表达量为0的基因剔除，不参与中位数计算，可以剔除特异性表达基因的影响；

3. 中位数算法可以进一步降低高差异表达基因的影响，而提高中等表达的基因的地位；

标准化因子的生物学意义

其实这个标准化因子算法就是选出一个有代表性的gene X（其实是每个样本一个代表性gene X），而这个gene X的reads for gene X/average for gene X比值就是标准化因子。

只不过选取gene X的时候，通过对数变换和中位数的方法，更多的参考了中等表达基因和管家基因的数据趋势，而剔除了特异性表达基因和高差异表达基因的影响。

相比较RPKM,FPKM,TPM标准化方法是除以总Read数，DESeq2标准化方法是除以一个有代表性基因的Read数，只不过这个Read数进行了变换（它除以了几何平均Read数， reads for gene X/average for gene X）。因为更能处理存在特异性表达基因和高差异表达基因的数据。

一个例子

我们按照这个DESeq2的标准化方法的思想，对图2中的数据进行一个简单的标准化，没有完全按照上述的7步法，只是体会这种标准化的意思（结果没有大的差异，但是算法并不完全正确）。

#Gene        Sample#1    Saple#2

A1BG        30          235
A1BG-AS1    24          188
A1CF        0           0
A2M            563         0
A2M-AS1        5           39
A2ML1        13          102

• 首先找到参照基因Gene X

  对于Sample 1#来说，A1BG-AS1和A2ML1是中位数基因（不计算A1CF、A2M），而样本 2#也是同样。

  为了简便，就用A2ML1基因了，其实结果是类似的。

原本应该根据ln(reads for gene X/average for gene X)的值计算中位数，这里直接根据原始Read值计算，会有一定的差异。

• 求解reads for gene X/average for gene X比值

  由于，

  average for gene A2ML1= (ln13+ln102)/2 = 3.12

  于是，

  reads for gene A2ML1/average for gene A2ML1= (13,102)/3.12 = (4.16,32.66)

  也就是说两个样本的标准化因子分别是4.16和32.66。

• 进行标准化变换

  将原始矩阵按照样本的不同除以各自的标准化因子，得下表：

  可以发现，除了A2M基因有表达差异外，其他基因表达无明显差异。

#Gene        Sample#1        Saple#2

A1BG        7.205869671     7.194095374
A1BG-AS1    5.764695737     5.755276299
A1CF        0               0
A2M            135.2301542     0
A2M-AS1        1.200978278     1.1939137
A2ML1        3.122543524     3.122543524

参考资料

StatQuest课程：https:///video-index/