前言——主要内容
这篇笔记是StatQuest系列视频教程的第59，60，61节。主要内容讲的是高通测序数据的差异基因分析，其中，第59节的内容是edgeR进行的文库均一化；第60节是DESeq2的文库均一化；第61节则是讲的是edgeR和DESeq2均一化的一些阈值选择。
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%百分位数
注：关于百分位数的一些内容，可以参见《StatQuest学习笔记06——分位数及其应用》这篇笔记。
对于每个样本，计算出校正后的数据的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。
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或者是其它的高通量测序技术通常会比较不同组织类型之间的测序数据差异，例如我们可能会比较肝脏与脾脏的差异。这个时候就可能会出现问题，因为在肝脏中，存在着某些特异性基因，这些基因只在肝脏中大量转录，而在脾脏中转录不活跃，这就是文库补偿造成的差异的一个案例，如下所示：
还有一种文库补偿造成差异的案例就是在同一种类型的细胞中会出现，例如针对某细胞，我们敲除掉这个细胞中的某个转录因子或者是某些具有调控功能的基因，那么野生型细胞和敲除型细胞的转录组就有可能出现差异，如下所示：
针对这个问题，我们看一个简单的案例，在这个案例中，我们有2个样本，这2个样本的总reads数是相同的，都是635个reads，我们假设它们的所有基因都是一样的，不过有一处不同，如下所示：
我们看一下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没有找到相应的中文译名，这里译为“不相关过滤”），也叫解决多重检测问题，如下所示：
有关假阳性
当我们做一个统计学检验时，就有一定的概率得到的是错误结论，关于这个问题，可以参见《StatQuest学习笔记11——p值详解》，简单来说就是，我们认为“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方法来解决，具体的可能参考《StatQuest学习笔记12——FDR及实现》这篇笔记，不过，即使是这样，我们还会遇到另外的一个问题，先看一个案例。
在下面的图形中，我们看到了两个独立的分布，其中红色曲线表示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个样本，此时要尤其注意异常值，如下所示：
