edgeR and DESeq2, part2-独立过滤

​ 过滤低read计数的基因，又名独立过滤(Independent Filtering)。

​ 每次我们做统计检验，都有可能得出错误的结论。简而言之，当我们说p值< 0.05显著，那么我们也可以说5%的情况下我们会报假阳性。

​ 当我们处理1或2个基因的差异表达时，这不是什么大问题，因为2次测试的5%是很小的，我们不太可能报告假阳性(false-positive)。然而，当我们检查基因组中的每一个基因(大约20000个)，看看哪些基因在癌细胞中被调控失调时，5%*20000=1000 false positive。好消息是，FDR和Benjamini-Hochberg方法弥补了这个问题,但是还存在问题，我们先看一个例子：

FDR过滤

假设我们有两个独立的分布，红色曲线表示小鼠品系X的体重，蓝色曲线表示小鼠品系Y的体重。如果我们测量3只X品系的老鼠，那些值很大可能接近红色分布中间。如果我们测量3只Y品系的老鼠，那些值很大可能接近蓝色分布中间。对这些体重进行t检验将导致p值< 0.05，我们将正确地得出测量结果来自不同的分布。

image-20210106110808050.png

然而，我们时不时会得到重叠的值。在这种情况下，p值将会很大。这是一个“假阴性”。

image-20210106110848581.png

现在，让电脑从这些分布中抽取1000个样本(3个X，3个Y)做1000个t检验。画1000个p值的直方图，我们得到949个真阳性(p值<0.05),5个假阳性(p值>0.05),我们做的每一个测试都应该是"真阳性"并且p值< 0.05，这是因为每一个测试都使用两个不同分布的样本。

image-20210106111606005.png

让我们通过添加没有区别的test来使它更真实。我们将从已经使用两个不同的分布计算的1000个p值开始。然后再加上1000个来自相同分布的样本值，这些p值应该是> 0.05，但偶尔(5%的时间)，我们会得到p值< 0.05。

image-20210106111840928.png

一共有993个p值小于0.05，949个真阳性来自第一组p值，44个假阳性来自第二组p值，因为仅仅只有4%的p值<0.05是假阳性，我们不需要使用FDR。但这只是因为我们编造了数据，如果是真的数据，我们不知道百分比，所以我们要用FDR。

image-20210106112434174.png

在使用FDR校正后，846个FDR调整后的p值仍然< 0.05。827个正阳性值剩下，为原949的89%，19个假阳性值剩下，占846个的2%。

image-20210106112928967.png

现在，让它更像RNA-Seq，p值的数目增加到6000。意味着，1000个样本来自不同的分布，5000个样本来自相同的分布，在大多数情况下来自相同样本应该会给出较大的p值。

image-20210106113246830.png

有1215个p值<0.05,949个p值是真阳性值，266个p值是假阳性值，22%的p值是假阳性值。

image-20210106113657633.png

FDR校正后，仅仅剩下256个p值<0.05,250个真阳性值剩下，占949的26%。6个假阳性值剩下，占256的2%。FDR在限制“显著结果”中的假阳性数量方面做得很好，但在保持真实阳性方面做得并不出色。

image-20210106114134742.png

现在让我们把样本量增加到11000，1000个样本来自不同的分布，10000个样本来自相同的分布，然后做11000次t检验，得到10000个p值。有1430个p值<0.05,949个p值是真阳性值，481个p值是假阳性值，34%的p值是假阳性值。FDR校正后，仅仅剩下56个p值<0.05,54个真阳性值剩下，占949的6%。2个假阳性值剩下，每当我们增加样本的检验数目时，通过FDR校正的真阳性的（p小于0.05）数目都会减少：

image-20210106115039256.png

下图中，绿线代表通过FDR校正后，真阳性p值<0.05的数量。橙线代表通过FDR校正后，假阳性p值<0.05的数量。这张图显示，尽管FDR可以控制假阳性的比率，随着检验数目的增多，真阳性的比例却在下降，这也表明Benjamini Hochberg方法还有改进的空间!

image-20210106115508646.png

edgeR和DESeq2过滤

​ 一般的想法是，具有超低read计数的基因不能提供有用信息，因此，它们可以从数据集中删除。换句话说，即使这些基因在生物学上是有趣的，如果在一种样本类型中只有1或2个转录本，而在另一种样本类型中只有3或4个转录本，就很难得到准确的read计数。

edgeR过滤

在做任何事情之前，edgeR建议去除所有的基因，除了那些在两个样品或更多样本中CPM>1的基因。CPM=Counts Per Million,它弥补了文库之间read深度的差异。

计算CPM:

计算公式：

image-20210106120653891.png

image-20210106120831864.png

现在我们有了所有样本中所有基因的CPM值，让我们去除所有的基因，除了那些在两个样品或更多样本中CPM>1的基因。

image-20210106121209831.png

edgR的方法很简单，但是你应该意识到测序深度会影响它。例如，如果一个样本有500万reads，CPM标准化因子=5000000/1000000=5,如果有5reads比对到一个基因，这个基因的CPM=5/5=1 CPM.如果一个样本有8000万reads，CPM标准化因子=80000000/1000000=80,此时1CPM=80reads。

image-20210106121555961.png

另一方面，有时你需要一个更大的CPM阈值，例如，你有50000reads比对到一个样本，标准化因子CPM=50000/1000000=0.05,如果你有一个read比对到一个基因上，它将变成1/0.05 = 20 CPM。即使这个基因是在生物学相关的水平上转录，因为你只能读到一个read，它仍然存在很大的噪音。

edgeR中CPM阈值

我们得到一个很好的阈值(cut off)？我们通过真实的数据集来说明

• 我从一位同事那里得到一个数据集，每个样本平均有2200万reads。(4个“野生型”及4个“敲除型”样本)
• 我在没有过滤单个基因的情况下对它进行了edgeR，这将生成原始的p值
• 使用不同的CPM阈值过滤掉基因，然后矫正p值

我们使用不同的CPM阈值过滤掉基因，然后矫正p值，绘制成图，其中x轴为最小CPM阈值，y轴表示经过FDR校正后，p值<0.05的基因数量。当x=0时，没有过滤掉任何基因，当x=1(推荐的阈值)时，此时基因数量为3400，因为有很多reads，建议的阈值太严格了。使用较低的阈值可以鉴别出大约200个更重要的基因。

image-20210106124901166.png

edgeR的主旨是：要小心，在计算p值后，尝试不同的CPM阈值(cutoff)。

DESeq和edgeR的区别

1. edgeR查看单个样本，并确保至少有2个CMP大于阈值。相比之下，DESeq2查看某个基因所有样本均一化reads的平均值，>阈值便保留。此时你可能会想:“酷!!我用DESeq2的方法，但是如果测量的基因有异常值呢?”

• DESeq2有一个异常值值检测方法(我们将在另一个StatQuest中讨论)，但它只在每个类别有两个以上的样本时才生效。

下图是分别使用DESeq和edgeR对同一批数据进行处理，它们都在同一区域达到峰值，这两种方法都会产生相似的阈值。现在让我们看看另一个不同之处。

2.另一个不同是改变了x轴，DESeq2绘制了分位数与显著基因数量的图，而不是最小CPM阈值。0%的基因低于阈值，20%的基因低于阈值，40%的基因低于阈值，分位数是有用的，因为正如我们所看到的，CPM依赖于测序深度，但无论如何，分位数总是分位数。无论库中有800万或8000万reads，10%的基因总是小于0.1分位数。

image-20210106131846588.png

我们可以选择分位数和最小CPM，

3.DESeq2对这些点拟合一条曲线，DESeq2然后在拟合曲线上找到最大位置。阈值是曲线上的最大位置，减去拟合曲线与原始值之间的标准差。换句话说，在峰值噪声范围内的第一个分位数是CPM阈值。如果没有原始值超过阈值，则不进行过滤。

image-20210106132824031.png

现在我们知道了edgeR和DESeq2是如何过滤基因的。

• edgeR是保留那些在2个或2个以上样本中，CPM大于最小阈值的基因
• DESeq2保留那些平均CPM大于最小CPM的基因，然后绘制显著基因与分位数的散点图，找到拟合曲线，再用最大值减去噪声，即是阈值

image-20210106133105442.png

建议

• 如果使用edgeR，在计算p值后计算CPM阈值
• 应用DEseq2的方法很容易找到最优CPM，以edgeR的基因选择标准
• 如果你这样做了，请确保你引用了这两个出版物!

• 如果你使用DESeq2，当每个分类只有2个样本时要小心异常值。

  
  
  image-20210106133521094.png