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

在转录组测序（RNA-Seq）中，基因的表达量是我们关注的重点。基因表达量的衡量指标有：RPKM、FPKM、TPM。

RPKM：是Reads Per Kilobase per Million mapped reads的缩写，代表每百万reads中来自于某基因每千碱基长度的reads数。RPKM是将map到基因的read数除以map到基因组上的所有read数(以million为单位)与RNA的长度(以KB为单位)。说实话，这个英文说明真的很费解，其实可以理解为“Reads Per Kilobase Per Million Reads”，即“每一百万条Reads中，对基因的每1000个Base而言，比对到该1000个base的Reads数”，计算公式。

FPKM：Fragments per Kilobase Million，FPKM意义与RPKM极为相近。二者区别仅在于，Fragment 与 Read。RPKM的诞生是针对早期的SE测序，FPKM则是在PE测序上对RPKM的校正。只要明确Reads 和 Fragments的区别，RPKM和FPKM的概念便易于区分。Reads即是指下机后fastq数据中的每一条Reads，Fragments则是指每一段用于测序的核酸片段，在SE中，一个Fragments只测一条Reads，所以，Reads数与Fragments数目相等；在PE中，一个Fragments测两端，会得到2条Reads，但由于后期质量或比对的过滤，有可能一个Fragments的2条Reads最后只有一条进入最后的表达量分析。总之，对某一对Reads而言，这2条Reads只能算一个Fragments，所以，Fragment的最终数目是Reads的1到2倍之间。

TPM：Transcripts Per Million，这个英文也很费解。先不纠结字面意思了，直接解释它的计算方法。TPM的计算分3步：

step1：根据基因/转录本长度校正count值；假设某基因count值为R1，则校正后count值为：R1/(L1/1000)； 【注: L1为该基因的长度】；

step2：计算total 校正后count值；即所有基因的校正后count值总和，Rtotal；

step3：计算TPM；TPM结果为：R1*1000*1000000/(L1*Rtoatl)。

RNA-seq是二代测序技术中用来表示基因表达量或丰度的方法。在衡量基因表达量时，若是单纯以map到的read数来计算基因的表达量，在统计上是不合理的。因为在随机抽样的情况下，序列较长的基因被抽到的机率本来就会比序列短的基因较高，如此一来，序列长的基因永远会被认为表达量较高，而错估基因真正的表现量，所以Ali Mortazavi等人在2008年提出以RPKM在估计基因的表现量。

edgeR与DESeq2这两种方法并不使用RPKM，FPKM，TPM等方法来进行均一化，edgeR与DESeq2在对文库进行均一化时要考虑两个方面的问题：

第一，测序深度（RPKM，FPKM，TPM方法也能做到）；

第二，文库补偿（library composition没有找到相应的中文译名，此处译为“文库补偿”），因为不同的样本含有不同的活跃基因（FPKM，FPKM和TPM做不到），如下所示：

edgeR均一化步骤

下面是edgeR进行文库均一化的步骤。

第一步：移除所有未转录的基因

我们先看下面的一批测序数据，在这批数据中，有3个样本，每个样本有5个基因（这个数据只是虚拟的，为了方便说明问题，实际测序中不可能只有这几个基因），如下所示：

其中我们可以发现，Gene5在这3个样本中的reads数都是0，因此我们要把Gene 5给去掉，如下所示：

第二步：选择参考样本

在这个步骤中，我们需要在一批样本中，挑选一个样本作为“参考样本”，随后，我们会使用这个样本来均一化其他的样本，我个人理解，这个参考样本其实就相当于qPCR中的参考基因，如下所示：

现在我们就遇到另外一个新问题了，什么是好的参考样本，什么是坏的参考样本，好坏的标准是什么？

我们先看一个案例，这个案例讲的是坏的参考样本，所有的样本如下所示，其中Sample #3是一个非常差的参考样本，Sample #3中所有基因的reads数的均一化都只是基于一个reads数，即Gene 3，因为其它的基因都是0，因此Sample #3有很大的噪音，如下所示：

为了避免找到这样极端的样本，edgeR会选择那些最“平均”的样本，如下所示：

我们现在看一下edgeR如何找到这个最“平均”的样本，我们再看一批数据，如下所示：

在寻找最“平均”的样本时，我们需要进行a，b，c，d这四步。

第a步：用总reads数校正每个样本

注：原文是Scale each sample by its total read counts，我这里使用“校正”来指scale。

计算出每个样本的所有基因的总reads数，如下图左图所示，然后使用每个样本中每个基因的reads数除以每个样本的总reads数，如下图右图所示：

第b步：计算75%百分位数

对于每个样本，计算出校正后的数据的75%百分位数的值，或者是小于75%百分位数的值，例如，对于样本1来说，它的75%百分位数是0.26，或者是小于0.26，如下所示：

对于样本2来说，它的75%百分位数是0.36，或者是小于0.36，如下所示：

对于样本3来说，它的75%百分位数是0.13，或者是小于0.13，如下所示：

现在把这3个样本的75%百分位数放在一起，如下所示：

第c步：计算平均75%百分位数

现在计算这3个样本的平均75%百分位数，加起来，除以3即可，如下所示：

第d步：找出最近接近于平均75%百分位数的样本

“参考样本”的标准就是它的75%百分位数最接近于平均75%百分位数，样本1，样本2和样本3的75%百分位数分别为0.26，0.36，0.13，它们与平均75%百分位数的差值分别为0.01，0.11，0.12，其中，最接近于0.26的样本是样本1，因此样本1就是“参考样本”，如下所示：

第三步：计算标准化因子

基本思路

在这一步骤中，我们需要选择一些基因集（genes，复数）来创建标准化因子（scaling factors），计算的过程就是针对“参考样本”，分别选择其余的样本的基因集来创建标准化因子（也就是说不同样本创建标准化因子的基因集不同）。我们在第二步中找到了参考样本，也就是Sample #1，现在分别针对Sample #1，来计算Sample #2和样本3的标准化因子，如下所示：

这是edgeR方法的一个局限，针对不同的样本，edgeR会使用不同的基因集来计算它们的标准化因子（原文是this is one of the ramifications of edgeR’s approach. Different samples use different genes to derive their scaling factors.不知道我理解的对不对，可以对比一下原文自己理解），如下所示：

现在我们以Sample #2为例，说明一下如何选择用于创建标准化因子的基因集，我们先来看一下用于选择的基因的不同类型，下图是一个XY二维坐标系，如下所示：

XY轴把这个坐标系分成了四个象限，这个很好理解，我们先看下图中的第1个点（蓝色的），这个点位于左上部分，它表示这个基因主要在参考样本（reference）中转录，如下所示：

我们再看第2个点，这个点位于最右的坐标轴上，根据坐标轴的说明，我们知道，这个基因主要在Sample #2中转录，如下所示：

我们再看第3个点（下图中红色圆圈标出的点），这个点位于正中间的坐标轴上，它表示这个基因的很多reads在参考样本和Sample #2中都存在，如下所示：

我们再看第4个点，它们于中间最下面的坐标轴上，根据坐标轴的说明，我们知道，在参考样本和Sample #2中，都很少含有这个基因的reads，如下所示：

这样，我们把所有的基因的reads数都放到这个坐标上，我们就可以发现，中间有很大一部分基因没有偏倚，没有偏倚的意思是说，在中间的这些基因，它们在参考样本和Sample #2中的reads数都差不多，区别不大，如下所示：

而edgeR就会选择中间区域（红色椭圆区域）里的这些基因，排除那些有偏倚的基因，如下所示：

计算过程

上面只是edgeR选择基因集的基本思路。现在我们看一下edgeR是如何计算标准化因子的。

第a步：过滤偏倚基因

计算所用的数据是已经均一化后的数据（也就是每个基因的reads数除以总的reads数），下图中表示的是全部的N个基因，基因虽然很多，但是这些基因均一化的方法还是与前面所述的4个基因是一样的，如下所示：

edgeR是通过倍数差异的log转换（log fold differences）来过滤偏倚基因（biased genes）的，它的公式是log fold differences = log2（Reference/sample 2)。从这个公式我们就可以看出来，如果Reference的某个基因的值相对于Sample #2比较高，那么这个log fold differences就是一个正值（图中指向3），如果Sample #2的某个基因的值相对于Reference比较高，那么这个log fold differences就是一个负值（图中指向-3），如下所示：

最终，我们会选择一个阈值，例如+/-6，超过这个范围就认为是极度偏倚的值，需要除去，如下所示：

我们现在以Gene 1为例来说明一下，它的计算公式如下所示：

由于参考样本的Gene1是0，因此这个最终的log fold differences数值为-Inf，它表示无穷小。现在我们看一下Gene2的计算，如下所示：

最终计算出所有的结果，如下所示：

此时，我们就要移除那些Inf的基因，也就是说，要移除那些在样本中（无论是同时在两个样本，还是任何一个样本中）reads数为0的基因，如下所示：

经过上面的计算，我们就得到了一个新的数据集，这个数据集是经过log fold differences转换后的数据集，此数据集用于排除偏倚基因。

第b步：计算几何均数

现在我们还需要另外一个数据集，用于观察哪些基因在两个样本是高转录和低转录的，如下所示：

计算基因的高转录和低转录时，首选要计算每个基因的几何均数（the geometric mean），几何均数很有用，因为它不太容易受到异常值的影响，如下所示：

这里要注明一下，严格意义来说，上面的这个公式只是经log转换后的数字的算术均数，由于经过了log转换，因此log转换后的算术均数其实上就是原始数据（未经log转换的数据）的几何均数，因为我们并不把log转换后的算术均数转换回正常的数值（也就是没有经过log转换后的数值），因此在这里，我们使用了log转换后的数值的算术均数，但它们的最终的效应是一样的，即异常值很少影响这些数据，下面以Gene 2为例说明一下这个计算过程：

Gene 2在参考样本中的数值为0.04，在Sample #2中的数值为0.05，它的几何均值其实应该是下面的这个样子：

现在对0.0447213595进行log2转换，如下所示：

我们再来看一下计算结果，如下所示：

从上面的结果我们可以发现，这个数值是一样的，经过计算，所有的结果如下所示：

现在我们移除那些infinite的值，也就是那些没有任何reads的基因（也就是计算结果中是-Inf或Inf的值），如下所示：

第c步：计算代表基因集

经过前面的计算，此时，我们就有了两张表，第一张表是log2(reference/Sample #2)的数据，它用于确定偏倚基因，另外一张表的数据是经log2转换后的均值数据，这批数据用于确定哪些基因是高转录的，哪些基因是低转录的，这两张表如下所示：

此时，我们将这两张表中的数据都按从小到大的顺序进行排列，如下所示：

在第一张表中，去掉前30%的数据，以及去掉后30%的数据，如下所示：

在第二张表中，去掉前5%的数据，以及去掉后5%的数据，如下所示：

用两张表中剩下的数据来计算标准化因子（取两张表基因的交集），如下所示：

不过在这个案例中，这两终表中剩余的基因并没有列出，都在省略号中，但是这个案例只是在讲edgeR的算法原理，并不涉及具体的数据，如下所示：

为了方便理解，我们就假定，下图就是最终于用计算Sample #2的标准化因子的基因集，如下所示：

同样的，还有计算Sample #3标准化因子的基因集，如下所示：

第四步：计算加权均数

在这一步骤中，edgeR会计算代表性基因集的log fold的加权平均数，在edgeR中，这个log fold的加权平均数称为为“weighted trimmed mean of the log2 ratio”，因为这些数据已经剔除掉（trim）了那些表达比较极端的基因（特别高的和特别低的），这样，计算结果就不会受到表达异常基因的影响，如下所示：

现在我们回到Sample #2的基因集上，这些基因是经过了第三步的筛选过的基因，只是我们假定的数据，就是我们上面说的省略号中的基因，只是用来计算的，不用管它的实际意义，它的log fold如下所示：

一旦我们选择了这些基因用于计算标准化因子，那么我们就仅需要计算出它们的log 2的加权均数即可，也就是说，比对到某个基因上的reads数越多，这个基因的权重就越大，这是因为对于那些低reads数的基因来的，log fold数值的变异程度比较大，如下所示：

为了说明reads数少的基因经过log fold转换后的变异程度比较大的问题，我们这里插入一个案例，下图是几个基因的reads数，如下所示：

从上面的数据中我们可以知道，Gene #1、Gene #2、Gene #3的reads数比较高，不过Gene #1比Gene #3的reads数还是多了4，如下所示：

我们再看Gene #4和Gene #6，虽然它们的reads数比较少，不过它们之间的reads数也是差了4，这个数值与Gene #1比Gene #3的差值一样，如下所示：

不过，我们看一下右侧的经过log2 fold转换后的数据，我们可以发现，Gene #1和Gene #3的log2 fold的值差别并不大，如下所示：

但是，我们看一下Gene #4和Gene #6，虽然它们的reads数相差是4，但是，它们的log2 fold的值相差却是1.6，这个值是非常大的值了，如下所示：

针对这种情况，我们就需要对那些有着高reads数的基因加上更大的权重，如下所示：

对于Sample #3的处理，也是如此，如下所示：

log2 fold的加权均数的计算公式如下所示：

简单来说就是log2 foldA的值乘以相应基因的reads数的和，然后除以所有的reads数之和。

第五步：将加权log2 fold值转换为真值

在这一步中，我们需要把前面过计算出来的加权平均值转换为真值（也就是log2转换前的数值）。

标准化因子的公式如下所示：

不过，此时的这个“标准化因子”并不是edgeR所使用的标准化因子，为了区分edgeR中的标准化因子，我暂时称它为“原始标准化因子”，如下所示：

按照上述方法，再计算出Sample #3的“原始标准化因子”，如下所示：

经过上面的计算，这个原始标准化因子还要进行进一步的计算，才是edgeR的标准化因子。

我们来看一个实际的案例，这个案例使用的一个RNA-Seq实验中的数据计算出来的原始标准化因子（就是按上面的公式计算出来的标准化因子），其中，WT2是参考样本，其余的样本根据WT2进行均一化，如下所示：

第六步：原始标准化因子的中心化

下图是我们计算出来的原始标准化因子，如下所示：

这些值在数轴上主要以0.95为中心进行排列，如下所示：

我们对这批数据进行“中心化”，其实就是用它们的每个值除以这4个数据的几何均数，几何均数就是如下所示：

最终计算的edgeR标准化因子（也就是中心化的原始标准化因子）如下所示：

虽然除以几何均数并不会改变最终的计算结果，但是Mark Robinson（edgeR包的开发者）说，经过中心化的处理就能赋予标准化因子一些数学上的美感。

此时，edgeR的标准化因子计算到此，接着把表达矩阵中的数值除以这个标准化因子即可，下一篇文章讲另外一种标准化文库的方法，即DESeq2。