跟着存档教程动手学RNAseq分析（五）：DESeq2基因水平差异表达分析

跟着存档教程动手学RNAseq分析（一）

跟着存档教程动手学RNAseq分析（二）

跟着存档教程动手学RNAseq分析（三）：使用DESeq2进行计数标准化

跟着存档教程动手学RNAseq分析（四）：使用DESeq2进行DE分析的QC方法

DESeq2差异表达分析

差异表达分析工作流的最后一步是将原始计数拟合到NB模型中，并对差异表达基因进行统计检验。在这一步中，我们主要想确定不同样本组的平均表达水平是否存在显著差异。

img

DESeq2论文[1]发表于2014年，但该软件包不断更新，可用于R通过Bioconductor。它建立在从DSS[2]和edgeR[3]方法中对离散度估计和广义线性模型的使用的良好思想上。

使用DESeq2进行差异表达分析涉及多个步骤，如下面的蓝色流程图所示。简单地说，DESeq2将对原始计数进行建模，使用标准化因子(大小因子)来解释库深度的差异。然后，它将估计基因的散度，并缩小这些估计，以产生更准确的分散估计，以建立计数模型。最后，DESeq2将拟合负二项模型，并使用Wald检验或似然比检验进行假设检验。

img

注：DESeq2由开发人员积极维护并不断更新。因此，重要的是要注意你正在使用的版本。最近，实施了一些相当大的变化，影响了产出。要了解更多关于对2014年原始论文中描述的方法所做的具体修改的细节，请参阅DESeq2文档中的这一节[4]。

关于DESeq2下的统计概念的其他细节在Rafael Irizarry的EdX课程“生命科学系列的数据分析”的材料[5]中很好地阐明。

运行DESeq2

在执行差异表达分析之前，通过QC期间的探索和/或先前的知识，了解数据中存在哪些变异来源是一个好主意。一旦你知道了变化的主要来源，你就可以在统计模型中分析或控制它们之前，通过将它们包括在你的设计公式中来消除它们。

例如，假设你有以下元数据：

img

如果你想检查不同治疗之间的表达差异，并且你知道主要的变化来源包括性别和年龄，那么你的设计公式将是：

design <- ~ sex + age + treatment

波浪线(~)应该始终处理你的因子，并告诉DESeq2使用公式对计数进行建模。注意，设计公式中包含的因素需要与元数据中的列名匹配。

复杂设计

DESeq2也允许分析复杂的设计。你可以通过在设计公式中指定来探索交互或“差异中的差异”。例如，如果你想探究性别对治疗效果的影响，可以在设计公式中指定如下：

design <- ~ sex + age + treatment + sex:treatment

由于相互作用项sex:treatment在公式中位于最后，DESeq2输出的结果将输出这一项的结果。

在DESeq2文档[6]中有对复杂设计的额外建议。此外，Limma文档[7]为创建更复杂的设计公式提供了额外的见解。

注意：需要帮助确定应该在元数据中显示哪些信息?我们有额外的材料[8]（TODO）强调批量rna序列规划的考虑。在开始实验之前，请先看一下这些材料，以便进行适当的实验设计。

MOV10差异表达分析

既然我们知道了如何向DESeq2指定模型，我们就可以对原始计数运行差异表达分析流程了。

要从原始计数数据中获得差异表达分析的结果，我们只需要运行2行代码!

首先，我们创建一个DESeqDataSet，就像我们在计数标准化一节中所做的那样，并指定包含原始计数的txi对象、元数据变量，并提供我们的设计公式：

## Create DESeq2Dataset object
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)

然后，要运行实际的差异表达分析，只需调用函数DESeq()。

## Run analysis
dds <- DESeq(dds)

通过将函数的结果重新分配回相同的变量名(dds)，我们可以填充我们的DESeqDataSet对象的槽。

img

从标准化计数数据到线性建模的一切都是通过使用一个单一的函数来完成的!这个函数将输出它执行的各个步骤的消息：

estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

注意：DESeq2中有一些单独的功能，允许我们以逐步的方式执行工作流中的每个步骤，而不是单个调用。我们演示了一个生成大小因子来创建标准化矩阵的示例。通过调用DESeq()，将为你运行每个步骤的单独函数。

DESeq2差异基因表达分析工作流程

通过以上两行代码，我们刚刚完成了DESeq2差异基因表达分析的工作流程。分析中的步骤输出如下：

img

我们将详细研究这些步骤中的每一个，以便更好地理解DESeq2是如何执行统计分析的，以及我们应该检查哪些指标来探索我们的分析质量。

第一步：估计size factor

差异表达分析的第一步是估计大小因子，这正是我们在标准化原始计数时所做的。

img

DESeq2将在执行差异表达分析时自动估计大小因子。但是，如果你已经使用estimateSizeFactors()生成了大小因子，就像我们前面所做的那样，那么DESeq2将使用这些值。

为了标准化计数数据，DESeq2使用前面在“计数标准化”一节中讨论的比值中值方法计算每个样本的大小因子。

MOV10 DE分析：检查尺寸因素

让我们快速看看每个样本的大小因子值：

> sizeFactors(dds)
Irrel_kd_1 Irrel_kd_2 Irrel_kd_3 Mov10_kd_2 Mov10_kd_3 Mov10_oe_1 Mov10_oe_2 Mov10_oe_3 
 1.1150371  0.9606366  0.7493552  1.5634128  0.9359082  1.2257749  1.1406863  0.6541689 

这些数字应该与最初运行函数estimateSizeFactors(dds)时生成的数字相同。看看每个样本的读段总数：

> colSums(counts(dds))
Irrel_kd_1 Irrel_kd_2 Irrel_kd_3 Mov10_kd_2 Mov10_kd_3 Mov10_oe_1 Mov10_oe_2 Mov10_oe_3 
  31160785   26504972   20498243   45300696   26745860   32062221   30025690   17285401 

这些数字与大小因素有什么关系？

现在来看看使用归一化后的总深度：

> colSums(counts(dds, normalized=T))
Irrel_kd_1 Irrel_kd_2 Irrel_kd_3 Mov10_kd_2 Mov10_kd_3 Mov10_oe_1 Mov10_oe_2 Mov10_oe_3 
  27945962   27591050   27354509   28975519   28577441   26156696   26322477   26423452 

如何将样本中的值与每个样本的总计数进行比较？

第二步：估计基因间散度

差异表达分析的下一步是估计基因间的离散。在我们进入细节之前，我们应该对DESeq2中离散指的是什么有一个很好的理解。

img

在RNA-seq计数数据中，我们知道：

1. 为了确定差异表达的基因，我们需要在给定组内（重复之间）差异的情况下，识别具有显著差异平均表达的基因。
2. 组内(重复之间)的变化需要考虑方差随平均表达增加的事实，如下图所示（每个黑点是一个基因）。

img

为了准确识别DE基因，DESeq2需要考虑方差和均值之间的关系。我们不希望所有的DE基因都是低计数的基因，因为低表达基因的方差更低。

它不使用方差作为数据变化的度量(因为方差与基因表达水平相关)，而是使用一种称为散度的变异度量，它解释了一个基因的方差和平均表达水平。离散度由Var = μ + α*μ^2计算，其中α表示离散度(Var =方差，μ =均值)，关系如下:

对于中等到高计数值的基因，散度的平方根将等于变异系数(Var / μ)。因此，0.01的离散度意味着在生物重复中，在平均预期值周围有10%的差异。具有相同均值的基因的离散估计只会根据它们的方差而不同。因此，散度估计反映了给定平均值（如下图中的黑点所示）在基因表达中的方差。

img

离散度和我们的模型有什么关系?

为了精确模拟测序计数，我们需要对每个基因生成组内变异（同一样本组重复之间的变异）的精确估计。由于每组只有几个（3-6）重复，每个基因的变异估计往往是不可靠的。

为了解决这个问题，DESeq2在基因之间共享信息，使用一种称为“收缩”的方法，根据该基因的平均表达水平生成更准确的变异估计。DESeq2假设具有相似表达水平的基因具有相似的离散度。

分别估计各基因的离散度：

DESeq2根据基因表达水平(平均重复计数)和方差估计每个基因的离散度。

第三步：对基因离散估计拟合曲线

工作流的下一步是拟合一条曲线到基因离散估计。数据拟合曲线背后的想法是，不同的基因将有不同规模的生物多样性，但在所有基因中，将有一个合理的离散估计的分布。

img

这条曲线在下面的图中显示为一条红线，它绘制了给定表达强度的基因的预期离散值的估计。每个黑点是一个具有相关的平均表达水平和离散的最大似然估计(MLE)的基因(步骤1)。

img

第四步：缩小基因间的离散度，使其接近曲线所预测的值

工作流的下一步是将基因的离散估计缩小到预期的离散值。

img

当样本量较小时，该曲线允许更准确地识别差异表达基因，每个基因的收缩强度取决于：

• 基因分散与曲线有多接近
• 样本量(样本量多=收缩量少)

这种收缩方法对于减少差异表达分析中的假阳性尤其重要。离散度估计低的基因向曲线收缩，输出更准确、更高的收缩值，用于模型拟合和差异表达检验。

略微高于曲线的离散估计也会向曲线收缩，以便更好地进行离散估计；然而，具有极高离散值的基因则不然。这是由于该基因可能不遵循建模假设，并且由于生物或技术原因比其他基因具有更高的可变性[1]。向曲线方向收缩值可能导致假阳性，因此这些值没有收缩。这些基因被下面的蓝色圆圈包围着。

img

这是一个很好的图，可以检查以确保你的数据很适合DESeq2模型。你期望你的数据通常在曲线上分散，随着平均表达水平的增加，这种分散会减少。如果你看到一片云或不同的形状，那么你可能想更多地探索你的数据，看看是否有污染(线粒体等)或异常值样本。请注意，对于任何低自由度的实验，在plotdispest()图中，你通过整个平均值范围获得了多少收缩。

令人担忧的分散图示例如下：

img

img

MOV10 DE分析:探索离散估计和评估模型拟合

让我们看看MOV10数据的离散度估计：

dds <- estimateDispersions(dds)
plotDispEsts(dds)

img

由于我们的样本量很小，对于许多基因，我们看到了相当大的收缩。你认为我们的数据适合这个模型吗?

第五步：对每个基因进行通用线性模型拟合

DESeq2工作流程的最后一步是为每个基因拟合负二项模型并进行差异表达检验。

img

如前所述，RNA-seq生成的计数数据表现出过度分散(方差>均值)，用于计数建模的统计分布需要考虑这种过度分散。DESeq2使用负二项分布来模拟RNA-seq计数，使用的公式如下：

img

建模是一种数学形式化的方法，可以在给定一组参数(即大小因子、离散度)的情况下，近似数据的行为。DESeq2将使用这个公式作为每个基因的模型，并将标准化计数数据拟合到其中。拟合模型后，估计各样本组的系数及其标准误差。这些系数是每个样本组的log2倍数变化的估计值。

负二项分布[9]

负二项分布（Negative binomial distribution）是统计学[10]上一种描述在一系列独立同分布的伯努利试验中，成功次数到达指定次数（记为r）时失败次数的离散概率分布[11]。“负二项分布”与“二项分布”的区别在于：“二项分布”是固定试验总次数N的独立试验中，成功次数k的分布；而“负二项分布”是所有到r次成功时即终止的独立试验中，失败次数k的分布。

img

img

假设检验

假设检验的第一步是为每个基因建立一个零假设。在我们的例子中，零假设是两个样本组之间没有差异表达(LFC == 0)。请注意，我们可以在不观察任何数据的情况下做到这一点，因为它是基于一个思想实验。其次，我们使用统计检验来确定根据观察数据，零假设是否为真。

Wald test

在DESeq2中，Wald检验是比较两组时默认使用的假设检验。Wald检验是对假设的检验，通常是对最大似然估计得到的参数进行检验。在我们的案例中，我们正在测试每个基因模型系数(LFC，log fold change)，它是通过使用最大似然估计利用散度等参数导出的。

DESeq2通过以下方式实现Wald检验：

• 用LFC除以标准误差，得到z统计量。
• 将z统计量与标准正态分布进行比较，并计算p值，报告随机选择出极端值至少为观测值的概率。
• 如果p值很小，我们拒绝零假设，并声明有证据反对零假设(即基因有差异表达)。

Wald Test and Z Test[12]

两者区别不大，但原文中说使用的是z统计量未必完全正确。

多重检验校正

如果我们直接使用来自Wald检验的p值，显著性临界值p < 0.05，这意味着有5%的可能性是假阳性。每个p值是单个检测(单个基因)的结果。我们测试的基因越多，假阳性率就越高。这就是多重检验问题。例如，如果我们测试20000个基因的差异表达，在p < 0.05的情况下，我们可能会偶然发现1000个基因。如果我们发现3000个基因有差异表达，大约三分之一的基因是假阳性。我们不想通过筛选我们的“重要”基因来确定哪些是真正的阳性基因。

DESeq2通过在检测前移除那些不太可能显著DE的基因，例如那些计数数量低和异常值样本(基因水平QC)，帮助减少被检测基因的数量。然而，我们仍然需要纠正多重检验，以减少假阳性的数量，有一些常见的方法：

• Bonferroni：调整后的p值通过以下方式计算：p值* m (m =测试总数)。这是一种非常保守的方法，错误否定的概率很高，因此通常不推荐使用。
• FDR/Benjamini-Hochberg: Benjamini和Hochberg(1995)定义了FDR的概念，并创建了一种算法，在给定一组独立的p值的情况下，将预期FDR控制在指定的水平以下。在DESeq2中，我们对控制FDR的BH方法进行了解释，我们将基因按p值排序，然后将每个排序后的p值乘以m/rank。
• q值/ Storey法：当该值显著时，可以达到的最小FDR。例如，如果基因X的q值为0.013，这意味着有1.3%的p值小于基因X的基因是假阳性。

那么FDR < 0.05是什么意思呢？通过将FDR截断值设置为< 0.05，我们表示，我们预期的差异表达基因的假阳性比例为5%。例如，如果您将500个基因称为差异表达，FDR截断值为0.05，那么预计其中25个是假阳性。

MOV10差异表达分析：控制组 VS 过表达组

我们有三个示例类，因此我们可以进行三种可能的成对比较：

1. Control vs. Mov10 overexpression
2. Control vs. Mov10 knockdown
3. Mov10 knockdown vs. Mov10 overexpression

我们只对上面的1和2感兴趣。利用设计公式，我们提供了~sampletype，表明这是我们感兴趣的主要因素。

创建对比contrasts

为了向DESeq2表明我们想要比较的两个组，我们可以使用对比。然后将对照品提供给DESeq2，使用Wald检验进行差异表达检验。DESeq2可以通过两种不同的方式提供对比：

1. 什么也不做。DESeq2将自动使用感兴趣条件的参考因子水平作为统计检验的基础。因子水平是根据级别的字母顺序选择的。
2. 在results()函数中，你可以指定感兴趣的比较和要比较的级别。最后给出的水平是进行比较的基础水平。语法如下所示：

# DO NOT RUN!
 contrast <- c("condition", "level_to_compare", "base_level")
 results(dds, contrast = contrast, alpha = alpha_threshold)

上面是示例，见下面实际操作。

构建results表格

为了构建结果表，我们将使用results()函数。为了告诉DESeq2我们想要比较哪些组，我们用contrast提供了我们想要做的对比。

## Define contrasts, extract results table, and shrink the log2 fold changes

contrast_oe <- c("sampletype", "MOV10_overexpression", "control")
res_tableOE <- results(dds, contrast=contrast_oe, alpha = 0.05)

上面我们提供了results()函数的最简单调用。看一下帮助手册，看看我们可以修改的其他参数：

?results

• 独立过滤：我们包括alpha参数并将其设置为0.05。这是用于优化独立滤波的显著性截止(默认设置为0.1)。如果我们想调整后的p值截止值(FDR)不是0.1(对于我们的最终重要基因列表)，alpha应该设置为该值。还有一个参数通过设置independentFiltering = F来关闭过滤。

什么是独立过滤?这是一个较低的平均阈值，是根据你的数据经验确定的，通过减少多重测试中被考虑的基因数量，可以增加显著基因的比例。

img

Image courtesy of slideshare presentation[13] from Joachim Jacob, 2014.

• Cooks cutoff：我们也可以通过cooksCutoff来过滤极端异常的基因。
• 多重校正：在DESeq2中，Wald检验得到的p值默认使用Benjamini和Hochberg方法进行多重校正。使用pAdjustMethod参数可以使用其他方法。

结果解释

结果表看起来非常像一个数据框，并且在许多方面可以像一个数据框架一样对待它(即在访问/子集数据时)。然而，重要的是要认识到它实际上存储在一个DESeqResults对象中。当我们开始将数据可视化时，这些信息将会很有帮助。

class(res_tableOE)

让我们浏览一下结果表中的一些列，以便更好地了解我们正在查看的内容。要提取关于每一列含义的信息，可以使用mcols():

mcols(res_tableOE, use.names=T)

• baseMean: mean of normalized counts for all samples
• log2FoldChange: log2 fold change
• lfcSE: standard error
• stat: Wald statistic
• pvalue: Wald test p-value
• padj: BH adjusted p-values

现在让我们看看结果中存储了哪些信息：

res_tableOE %>% data.frame() %>% View()
log2 fold change (MAP): sampletype MOV10_overexpression vs control 
Wald test p-value: sampletype MOV10_overexpression vs control 
DataFrame with 57914 rows and 6 columns
                 baseMean log2FoldChange lfcSE  stat  pvalue  padj
                <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
ENSG00000000003  3.53E+03 -0.427190489 0.0755347 -5.65604739 1.55E-08 4.47E-07
ENSG00000000005  2.62E+01 0.016159765 0.23735203 0.06584098 9.48E-01 9.74E-01
ENSG00000000419  1.48E+03 0.362663551 0.10761742 3.36995355 7.52E-04 4.91E-03
ENSG00000000457  5.19E+02 0.219135591 0.09768842 2.24476439 2.48E-02 8.21E-02
ENSG00000000460  1.16E+03 -0.261603812 0.07912962 -3.30661411 9.44E-04 5.92E-03
...   ...  ...  ...

对比中名称的顺序决定了所报告的foldchange变化的方向。对比到第二个元素中提供的名称是用作baseline的水平。例如，如果我们观察到-2的log2倍的变化，这意味着相对于对照组，Mov10_oe中的基因表达更低。然而，这些估计并不能解释我们在低读计数情况下观察到的巨大离散。为了避免这种情况，需要调整模型计算的log2 fold change。

虽然知道所提供的fold changes是重要的，但最终应使用p调整值来确定重要的基因。显着基因可以输出为可视化和/或功能分析。

注意：p值设置为NA

1. 如果在一行中，所有样本计数为零，则baseMean列将为零，log2倍的变化估计值、p-value和调整后的p-value都将设置为NA。
2. 如果一行包含一个极端的计数异常值的样本，那么p-value和调整后的p-value将被设置为NA。这些异常值是通过库克距离来检测的。
3. 如果一行被自动独立过滤，由于其归一化计数平均值较低，则只有调整后的p值将被设置为NA。

搜索log2 fold change (LFC)

为了生成更准确的log2 fold change估计，DESeq2允许对基因信息较低时LFC估计缩窄至零，其中可能包括：

• 低计数
• 高散度值

与离散度估计的收缩一样，LFC收缩利用来自所有基因的信息来产生更准确的估计。具体来说，所有基因的LFC估计的分布(作为一种先验)被用来缩小信息少或高度分散的基因的LFC估计，向更可能(更低)的LFC估计。

img

Illustration taken from the DESeq2 paper[14].

例如，在上图中，绿色基因和紫色基因在两个样本组(C57BL/6J和DBA/2J)中有相同的平均值，但绿色基因变异少，紫色基因变异水平高。对于低变异的绿色基因，未收缩的LFC估计(绿色实线顶点)与收缩的LFC估计(绿色虚线顶点)非常相似，但对于紫色基因，由于高度分散，LFC估计差异较大。因此，尽管两个基因可以有相似的标准化计数值，但它们的LFC收缩程度可能不同。请注意，LFC估计向先验方向(黑色实线)缩小了。

在DESeq2的最新版本中，默认情况下不会执行LFC估算的缩减。这意味着log2折叠变化将与通过以下方法计算的结果相同：

log2 (normalized_counts_group1 / normalized_counts_group2)

为了生成缩窄的log2倍更改估计，你必须使用函数lfcShrink()在结果对象(我们将在下面创建)上运行额外的步骤。

## Save the unshrunken results to compare
res_tableOE_unshrunken <- res_tableOE

# Apply fold change shrinkage
res_tableOE <- lfcShrink(dds, contrast=contrast_oe, res=res_tableOE)

注：将log2 fold change缩小并不会改变被鉴定为显著差异表达的基因的总数。fold change的收缩是为了帮助下游评价结果。例如，如果你想要根据你的重要基因的fold change来进行进一步的评估，你可能想要使用缩小后的值。此外，对于像GSEA这样的功能分析工具，它需要fold change值作为输入，而你希望提供缩小后的值。

MA图

一个对探索我们的结果有用的图是MA图。MA图显示了所有基因测试的归一化计数的平均值与log2 fold change。显著DE的基因被着色以便于识别。这也是一个很好的方式来说明LFC收缩的影响。DESeq2包提供了一个生成MA图的简单函数。

让我们先绘制没有收缩的结果：

plotMA(res_tableOE_unshrunken, ylim=c(-2,2))

img

然后绘制收缩后的结果：

plotMA(res_tableOE, ylim=c(-2,2))

img

除了上面描述的比较之外，这个图让我们可以评估fold change的幅度，以及它们相对于平均值的分布情况。一般来说，我们希望在整个表达水平范围内看到重要的基因（而不是偏向表达量高或低的基因）。

MOV10差异表达分析：Control vs Knockdown

现在我们有了过表达的结果，让我们对对照和敲除样本做同样的操作。使用results()中的contrast提取结果表并将其存储到名为res_tableKD的变量中。

## Define contrasts, extract results table and shrink log2 fold changes
contrast_kd <-  c("sampletype", "MOV10_knockdown", "control")

res_tableKD <- results(dds, contrast=contrast_kd, alpha = 0.05)

res_tableKD <- lfcShrink(dds, contrast=contrast_kd, res=res_tableKD)

快速浏览一下结果表，其中包含我们感兴趣的Control-Knockdown比较的Wald检验统计数据，并确保其格式与我们在OE中观察到的类似。

汇总结果

为了对结果表进行汇总，DESeq2中的一个方便的函数是summary()。令人困惑的是，它与用于检查数据框的函数同名。当以DESeq结果表作为输入调用该函数时，将使用alpha阈值:FDR < 0.05(即使输出显示p-value < 0.05，仍然使用padj/FDR)汇总结果。让我们从OE vs控制结果开始：

## Summarize results
summary(res_tableOE, alpha = 0.05)

除了在默认阈值上上调和下调的基因数量，该函数还报告被检测的基因数量(总读计数非零的基因)，以及由于平均计数低而没有包括在多次检测校正中的基因数量。

提取显著差异表达基因

让我们首先创建包含阈值条件的变量：

### Set thresholds
padj.cutoff <- 0.05

我们可以使用filter()函数轻松地提取结果表子集，只包含那些重要的结果，但首先我们要将结果表转换为tibble：

res_tableOE_tb <- res_tableOE %>%
  data.frame() %>%
  rownames_to_column(var="gene") %>% 
  as_tibble()

现在我们可以用我们预先设定的阈值来划分这个表，只保留重要的基因：

sigOE <- res_tableOE_tb %>%
        filter(padj < padj.cutoff)
sigOE

使用与上面相同的p调整阈值(padj。cutoff < 0.05)，子集res_tableKD报告与对照组相比，Mov10_knockdown中上调和下调的基因数量。

res_tableKD_tb <- res_tableKD %>%
  data.frame() %>%
  rownames_to_column(var="gene") %>% 
  as_tibble()
  
sigKD <- res_tableKD_tb %>%
        filter(padj < padj.cutoff)

与对照相比，敲低基因中有多少基因有差异表达?

sigKD

既然我们已经提取了重要的结果，我们就可以进行可视化了!

添加一个fold change阈值

有大量的重要基因列表，很难提取有意义的生物相关性。为了帮助增加严格性，还可以添加fold change阈值。

例如，我们可以创建一个新的阈值lfc。切断并将其设置为0.58(请记住，我们使用的是log2倍的更改，因此这可以转换为实际的1.5倍更改)。

添加fold change阈值的另一种方法：

results()函数有一个使用lfcThrehsold参数添加fold change阈值的选项。这种方法是统计学驱动，当你想要一组基于特定fold change阈值的更自信的基因时，推荐使用这种方法。它实际上对所期望的阈值执行一个统计检验，即对大于指定绝对值的log2倍的变化执行一个双尾检验。用户可以使用假设来改变备择假设，也可以进行两次单侧检验。这是一种更保守的方法，所以期望检索到更小的一组基因!

参考资料

[1]

DESeq2论文: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8

[2]

DSS: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4005660/

[3]

edgeR: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2796818/

[4]

DESeq2文档中的这一节: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#methods-changes-since-the-2014-deseq2-paper

[5]

“生命科学系列的数据分析”的材料: https://rafalab.github.io/pages/harvardx.html

[6]

DESeq2文档: https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions

[7]

Limma文档: http://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf

[8]

额外的材料: https://hbctraining.github.io/Intro-to-rnaseq-hpc-salmon/lessons/experimental_planning_considerations.html

[9]

负二项分布: https://zh.m.wikipedia.org/zh-hans/负二项分布

[10]

统计学: https://zh.m.wikipedia.org/wiki/統計學

[11]

离散概率分布: https://zh.m.wikipedia.org/wiki/离散概率分布

[12]

Wald Test and Z Test: https://stats.stackexchange.com/questions/152630/wald-test-and-z-test

[13]

slideshare presentation: https://www.slideshare.net/joachimjacob/5rna-seqpart5detecting-differentialexpression

[14]

DESeq2 paper: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8