edgeR文库均一化详细过程

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

        以下面的虚拟测序数据为例，在这批数据中，有3个样本，每个样本有5个基因，如下所示：

可以看到，Gene5在这3个样本中的reads数都是0，因此要把Gene 5去掉，如下所示：

第二步：选择参考样本
        我们需要在这组样本中，挑选出一个样本作为“参考样本”，然后使用这个样本来均一化其他的样本。如果在一个样本中只有一个基因的read数大于0，其他基因的read全是0，选择这个样本作为“参考样本”，会导致所有基因的reads数的均一化都只是基于一个reads数，这会引入很大的噪音。为了避免出现这种极端情况，edgeR会选择最“平均”的样本作为“参考样本”。

挑选参考样本步骤：

（1）用总reads数校正每个样本

       分别计算出每个样本中计算75%百分位数所有基因的总reads数，如下图左图所示，然后把每个样本中每个基因的reads数除以相应样本的总reads数，如下图右图所示：

（2）计算75%百分位数

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

       样本2和样本3的75%百分位数计算方法跟样本1相同。把这3个样本的75%百分位数放在一起，如下所示：

 （3）计算平均75%百分位数

        把3个样本的75%百分位数加起来除以3，即：(0.26+0.36+0.13)/3 = 0.25。

 （4）找出最近接近于平均75%百分位数的样本

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

第三步：计算标准化因子

基本思路

      在这一步骤中，我们需要选择一些基因集来创建标准化因子（scaling factors），计算的过程就是针对“参考样本”，分别选择其余样本的基因集来创建标准化因子（针对不同的样本，edgeR会使用不同的基因集来计算它们的标准化因子）。在选择基因集的时候，edgeR会选择那部分没有偏倚的基因，没有偏倚的意思是说，这些基因在参考样本和Sample #2中的reads数都差不多，区别不大。

计算过程

（1）过滤偏倚基因

       计算所用的数据是已经把样本中每个基因的reads数除以对应样本总的reads数以后的数据。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，超过这个范围就认为是极度偏倚的值，需要除去。

         现在以Sample #2为例来展示一下计算过程，Gene1：

 Gene2:

       计算出样本2所有的结果，并且移除那些Inf的基因，也就是说，移除那些在样本中（无论是同时在参考样本和样本2中，还是任何一个样本中）reads数为0的基因，如下所示：

       经过上面的计算，我们就得到了一个新的数据集，这个数据集是经过log fold differences转换后的数据集，此数据集用于排除偏倚基因（超过指定阈值范围的基因就认为是偏倚基因，删除）。

（2）计算几何均数

       现在我们还需要另外一个数据集，用于观察哪些基因在两个样本是高转录和低转录的，用到的数据跟步骤（1）相同，是已经把样本中每个基因的reads数除以对应样本总的reads数以后的数据。计算基因的高转录和低转录时，首先要计算每个基因的几何平均数（the geometric mean），几何平均数很有用，因为它不太容易受到异常值的影响，如下所示：

 

        这里要说明一下，严格意义来说，上面的这个公式只是经log2转换后的数字的算术平均数，但是数据经过了log转换，根据log的特殊运算法则可知，此时log转换后的算术平均数实际上就是原始数据（未经log转换的数据）的几何平均数，取log2。

        通过上述方式计算出Sample #2中所有基因与参考样本基因的几何平均数，如下图所示：

        接下来移除那些-Inf或Inf的值，也就是那些没有任何reads的基因，如下所示：

 （3）计算代表基因集

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

       把这两张表中的数据按从小到大的顺序排列，然后去掉第一张表的前30%和后30%的数据，去掉第二张表的前5%和后5%的数据。取两张表基因的交集，用于后续计算样本2的标准化因子，如下图：

经过这一步，我们已经得到了样本2用于计算标准化因子的代表基因。

（4）计算加权平均数 

       在这一步骤中，edgeR会计算代表性基因集的log fold的加权平均数，在edgeR中，这个log fold的加权平均数称为为“weighted trimmed mean of the log2 ratio”，因为这些数据已经剔除掉（trim）了那些表达比较极端的基因（特别高的和特别低的），这样，计算结果就不会受到表达异常基因的影响。通过前面3步，我们得到了用于计算标准化因子的代表基因集，此时，仅需计算出它们的log2(Reference/sample 2)的加权平均数即可。log2 fold的加权平均数的计算公式如下所示：

        简单来说就是log2 fold的值乘以相应基因的reads数的和，然后除以样本中所有基因的reads数之和。注意：这里所说的log2 fold值是两张表基因交集的log2(Reference/sample 2)。

 （5）将加权log2 fold值转换为真值

        在这一步中，我们需要把前面过计算出来的加权平均值转换为真值（也就是log2转换前的数值），以2为底数，取log2 fold次方即可，此时得到了Sample #2最终的原始标准化因子。这里的原始标准化因子还不是egdeR最终使用的标准化因子。

（6）原始标准化因子的中心化

       重复以上5个步骤，我们可以得到所有样本的原始标准化因子，最后要对得到的原始标准化因子进行“中心化”。首先计算原始标准化因子的几何平均数，然后把每个原始标准化因子除以几何平均数。此时，得到了egdeR使用的最终的标准化因子。

第四步：样本文库校正

        综合上述三步，已经得到了edgeR的标准化因子，把原始表达矩阵中的数值按列对应除以标准化因子就得到了文库校正后的表达矩阵。