2021-09-16差异表达基因时的Log2FC和FDR值的含义？

Posted 2023-05-02

参考技术A 转录组分析 差异表达 基因时，结果中会出现 Log2FC 和 FDR 值，这两个是什么意思呢？

log2FC 中的FC即 fold change，表示两样品（组）间表达量的比值，对其取以2为底的对数之后即为log2FC。一般默认取log2FC绝对值大于1为差异基因的筛选标准；

FDR 即False Discovery Rate，错误发现率，是通过对差异显著性p值（p-value）进行校正得到的。由于转录组测序的 差异表达 分析是对大量的基因表达值进行独立的统计假设检验，会存在假阳性问题，因此在进行差异表达分析过程中，采用了公认的Benjamini-Hochberg校正方法对原有假设检验得到的显著性p值（p-value）进行校正，并最终采用FDR作为差异表达基因筛选的关键指标。一般取FDR<0.01或者0.05作为默认标准。

这两个指标的选取一般是按照经验值去筛选的，并非完全不可以调整。在实验差异基因数目过低或者过高，可以对指标进行微调。

实际上经常看到的差异表达火山图（如下图）里的几条虚线就是这两个指标的体现。

为什么要用FDR

在转录组分析中，如何确定某个转录本在不同的样品中表达量是否有差异是分析的核心内容之一。一般来说，我们认为，不同样品中，表达量差异在两倍以上的转录本，是具有表达差异的转录本。为了判断两个样品之间的表达量差异究竟是由于各种误差导致的还是本质差异，我们需要根据所有基因在这两个样本中的表达量数据进行假设检验。常用的假设检验方法有t-检验、卡方检验等。很多刚接触转录组分析的人可能会有这样一个疑问，一个转录本是不是 差异表达 ，做完假设检验看P-value不就可以了么？为什么会有FDR这样一个新的概念出现？这是因为转录组分析并不是针对一个或几个转录本进行分析，转录组分析的是一个样品中所转录表达的所有转录本。所以，一个样品当中有多少转录本，就需要对多少转录本进行假设检验。这会导致一个很严重的问题，在单次假设检验中较低的假阳性比例会累积到一个非常惊人的程度。举个不太严谨的例子。

假设现在有这样一个项目：

● 包含两个样品，共得到10000条转录本的表达量数据，

● 其中有100条转录本的表达量在两个样品中是有差异的。

● 针对单个基因的 差异表达 分析有1%的假阳性。

由于存在1%假阳性的结果，在我们分析完这10000个基因后，我们会得到100个假阳性导致的错误结果，加上100条真实存在的结果，共计200个结果。在这个例子中，一次分析得到的200个差异表达基因中，有50%都是假阳性导致的错误结果，这显然是不可接受的。为了解决这个问题，FDR这个概念被引入，以控制最终得到的分析结果中假阳性的比例。

如何计算FDR

FDR的计算是根据假设检验的P-value进行校正而得到的。一般来说，FDR的计算采用Benjamini-Hochberg方法（简称BH法），计算方法如下：

1. 将所有P-value升序排列.P-value记为P，P-value的序号记为i，P-value的总数记为m

2. FDR(i)=P(i)*m/i

3. 根据i的取值从大到小，依次执行FDR(i)=minFDR(i),FDR(i+1)

注：实际上，BH法的原始算法是找到一个最大的i，满足P≤i/m*FDR阈值，此时，所有小于i的数据就都可以认为是显著的。在实践中，为了能够在比较方便的用不同的FDR阈值对数据进行分析，采用了步骤3里的方法。这个方法可以保证，不论FDR阈值选择多少，都可以直接根据FDR的数值来直接找到所有显著的数据。

下面我们以一个包含10个数据的例子来看一下FDR计算的过程

在这个例子中，第一列是原始的P-value，第二列是排序后的序号，第三列是根据P-value校正得到的初始FDR，第四列是最终用于筛选数据的FDR数值。如果我们设定FDR<0.05，那么绿色高亮的两个数据就是最终分析认为显著的数据。

FDR的阈值选择在转录组分析中是非常重要的一个环节，常用的阈值包括0.01、0.05、0.1等。实践中也可以根据实际的需要来灵活选择。例如，在做真菌或者原核生物的转录组分析时，由于这些物种转录本数量较少，假阳性累积的程度较低，所以可以适当将FDR阈值设置的较高一些，这样可以获得较多的差异表达结果，有利于后续的分析。

简单的转录组差异基因表达分析 -- DESeq2

参考技术A

经典的转录组差异分析通常会使用到三个工具 limma/voom , edgeR 和 DESeq2 。今天我们就通过一个小规模的转录组测序数据来演示 DESeq2 的简单流程。

对于 DESeq2 的分析流程而言，我们需要输入的数据包括：

​ 下面就以 mobData 中的数据为例简单介绍 DESeq2 的分析流程

​ 由于 mobData 中的行名没有提供基因的ID，我们也不是为了探究生物学问题，就以 mobData 的行数作为其ID

DESeqDataSet 是 DESeq2 流程中储存read counts和中间统计分析数据的对象，之后的分析都建立在该对象之上进行。

​ 在进行差异分析之前，需要对样本数据的表达矩阵进行预处理，包括：

​ 通过PCA结果来看各组样本分组情况还是不错的，但hclust的聚类结果反映的分组就略微有点混杂了，可能要聚类计算的距离函数选用不当有关。

使用 DESeq() 函数进行差异分析时，该函数干了以下三件事：

​ counts() 可以提取 DESeq object 中的表达矩阵，而 results() 可以提取差异分析的结果，其中包括了：

样本间的均值, log2 fold changes, standard errors, test statistics, p-values and adjusted p-values.

​ 使用 results() 函数时需要指明进行比较的样本，这里用 contrast=c("group_list","MM","WW") 提取 MM 组和 WW 组进行差异分析的结果。如果想要比较 WM 组和 WW 组，只要改变 contrast=c("group_list","WM","WW") 即可。

​ 检查结果中是否包含 NA 值

​ 这里 padj 中有1142个 NA 值是因为使用 results() 提取差异分析结果时，大于 alpha 值（这里是0.1）的矫正后p-value都会被当做是 NA 。因此，我们将这些 padj 值都设为 1

排序后以 log2FoldChange 绝对值大于1， padj 小于0.05为条件筛选显著的差异表达基因

​ 至此，便筛选出了 217个 在 MM 组和 WW 组之间的显著差异表达基因。至于后续的可视化分析则是因课题而异了，等以后有空了再补坑吧！

​ 有同学可能注意到，虽然我们的样本有多个组，但在差异分析时进行的还是pairwise的分析，为什么我们不可以三个组一起分析呢？学过ANOVA的同学应该都知道，ANOVA就是可以应对这种多组差异分析的情况。但要注意的是，ANOVA只可以告诉我们对于某个基因在这三组中是否存在差异，想要找出是哪一组有其他组别有差异还是需要进行pairwise t-test之类的分析。所以，在这里我们两组两组地进行分析正是出于这个考虑，并且有更方便我们解释差异分析的结果，说明在A组基因的表达量相对于B组的是上调还是下调。另外，本文的差异分析还是处于单因子水平（只有一个变量），至于多因子的差异分析以后研究透了再和大家进行分享。

完。