文章目录

学习目标
探索结果（Wald test）
- 指定的对比（Specifying contrasts）
- 我选择什么作为base level有关系吗?
- 结果表
- P值
- - 基因水平过滤（Gene-level filtering）
- 倍数变化（Fold change）
- - 更准确的LFC估计
- MA图

学习目标

讨论为两两比较(Wald检验)生成结果表所需的步骤
总结不同水平的基因过滤
解释log倍数变化收缩（log fold change shrinkage）

探索结果（Wald test）

默认情况下，DESeq2使用Wald检验来识别两个样本类之间差异表达的基因。考虑到设计公式中使用的因子，以及存在的因子水平的数量，我们可以提取一些不同比较的结果。在这里，我们将介绍如何从dds对象获取结果，并就如何解释这些结果提供一些解释。

注意:Wald检验也可以用于连续变量。如果在设计公式中提供的感兴趣的变量是连续值，那么报告的log2FoldChange 是该变量的每单位变化。

指定的对比（Specifying contrasts）

在我们的数据集中有三个样本类，所以我们可以进行三种可能的两两比较:

Control vs. Mov10 overexpression
Control vs. Mov10 knockdown
Mov10 knockdown vs. Mov10 overexpression

我们只对上面的第一条和第二条感兴趣。当我们创建dds对象时，我们提供了~ sampletype作为设计公式，表明sampletype是我们感兴趣的主要因素。
为了指出我们想要比较的两个样本类别，我们需要指定对比（contrasts）。对比用作DESeq2 results()函数的输入，以提取所需的结果。
对比可以用两种不同的方式来表示(第一种更常用):

对比可以作为一个字符向量，包含三个元素:设计公式中(感兴趣的)因素的名称，两个要比较的因素层次的名称。最后给出的因素水平是进行比较的基础水平。语法如下:

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

对比可以以2个字符串向量的列表：感兴趣基因倍数变化水平的名称，基线水平倍数变化名称。这些名称应该与resultsNames(object)的元素相同。该方法可用于组合交互项和主效应。

# DO NOT RUN!resultsNames(dds) # to see what names to usecontrast <- list(resultsNames(dds)[1], resultsNames(dds)[2])results(dds, contrast = contrast)

或者，如果你只有两个因素水平，你可以什么也不做，而不用担心指定对比(即results(dds))。在本例中，DESeq2将根据level的字母顺序选择你的基础因子水平。
首先，我们要计算MOV10过表达样本和对照样本之间的表达变化。因此，我们将使用第一种方法来进行对比，并创建一个字符向量:

## Define contrasts for MOV10 overexpression
contrast_oe <- c("sampletype", "MOV10_overexpression", "control")

我选择什么作为base level有关系吗?

是的，这很重要。决定哪个level是base level将决定如何解释所报告的fold change。例如，如果我们观察到-2的log2倍的变化这意味着基因表达在相关因素水平上比基础水平更低。因此，如果把它留给DESeq2来决定对比，一定要检查字母顺序是否与你预期的fold change方向一致。

结果表

现在已经创建了对比（contrast），可以使用它作为results()函数的输入。让我们快速看一下这个功能的帮助手册:

?results

你将看到，我们可以选择提供大量的参数，并根据需要调整默认值。在我们学习这节课的过程中，我们会不断回到帮助文档来讨论一些值得了解的参数。

## Extract results for MOV10 overexpression vs control
res_tableOE <- results(dds, contrast=contrast_oe, alpha = 0.05)

注意:对于我们的分析，除了contrast参数外，我们还将为alpha参数提供一个0.05的值。当我们谈到基因水平过滤（gene-level filtering）时，将更详细地描述这一点。

返回给我们的结果表是一个DESeqResults对象，它是DataFrame的一个简单子类。在许多方面，它可以被视为一个数据框(即在访问/提取子集数据时)，但是重要的是要认识到在可视化等下游步骤中存在差异。

# Check what type of object is returned
class(res_tableOE)[1] "DESeqResults"
attr(,"package")
[1] "DESeq2"

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

# What is stored in results?
res_tableOE %>%
data.frame() %>%
View()

log2 fold change (MLE): sampletype MOV10_overexpression vs control
Wald test p-value: sampletype MOV10 overexpression vs control
DataFrame with 57761 rows and 6 columnsbaseMean log2FoldChange     lfcSE       stat      pvalue        padj<numeric>      <numeric> <numeric>  <numeric>   <numeric>   <numeric>
ENSG00000000003  3525.8835      -0.438245 0.0774607 -5.6576466 1.53463e-08 4.25097e-07
ENSG00000000005    26.2489       0.029208 0.4411295  0.0662119 9.47209e-01 9.72687e-01
ENSG00000000419  1478.2512       0.383635 0.1137609  3.3722920 7.45454e-04 4.67392e-03
ENSG00000000457   518.4220       0.228971 0.1023313  2.2375412 2.52510e-02 8.02348e-02
ENSG00000000460  1159.7761      -0.269138 0.0814993 -3.3023359 9.58832e-04 5.76086e-03
...                    ...            ...       ...        ...         ...         ...
ENSG00000285889    1.82171       -4.68144 3.9266061   -1.19224  0.23316908          NA
ENSG00000285950    7.58089       -1.01978 1.0715583   -0.95168  0.34125953          NA
ENSG00000285976 4676.24904        0.19364 0.0656673    2.94880  0.00319006   0.0157503
ENSG00000285978    2.25697        4.13612 2.0706212    1.99753  0.04576795          NA
ENSG00000285980    0.00000             NA        NA         NA          NA          NA

我们为每个基因(行)报告了六列信息。我们可以使用mcols()函数来提取每个列中存储的值所代表的信息:

# Get information on each column in results
mcols(res_tableOE, use.names=T)

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

P值

p值是一个概率值，用来确定是否有证据拒绝零假设。较小的p值意味着有更有力的证据支持备择假设。 然而，因为我们是对每个个体基因进行测试，我们需要纠正这些进行多次检验的p值。
结果表中的padj列表示多次检验调整后的p值，是结果中最重要的一列。通常，阈值(如padj < 0.05)是识别显著基因的良好起点。DESeq2中多次检验校正的默认方法是Benjamini-Hochberg错误发现率(FDR)。还有其他可用的更正方法，可以通过向results()函数添加pAdjustMethod参数来更改。

p.adjust.methods
# c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
#   "fdr", "none")

基因水平过滤（Gene-level filtering）

让我们仔细看看我们的结果表。当我们滚动它时，您将注意到对于所选的基因，在pvalue和padj列中有NA值。这意味着什么?

缺失值表示作为DESeq()函数的一部分经过过滤的基因。在进行差异表达分析之前，剔除那些很少或没有机会被检测到差异表达的基因是有益的。这将增加检验出差异表达基因的效能。DESeq2不会从原始计数矩阵中直接移除任何基因，所以所有的基因都将出现在你的结果表中。DESeq2剔除的基因符合以下三个过滤标准之一:

所有样本中计数为0的基因
如果在一行中，所有的样本都是零计数没有表达信息，则这些基因不被检验。

# Filter genes by zero expression
res_tableOE[which(res_tableOE$baseMean == 0),] %>%
data.frame() %>%
View()

这些基因的baseMean列为零，log2倍数变化估计值、p值和调整后p值都将设置为NA。

基因计数为极端离群值
DESeq()函数为每个基因和每个样本计算一个被称为库克距离(Cook 's distance)的异常值诊断测试。库克距离是衡量单个样本对基因拟合系数的影响程度，库克距离的较大值是用来表示离群值。包含超过阈值的库克距离的基因被标记，但是至少需要3次重复来标记，因为很难判断哪个样本可能是只有2次重复的离群值。我们可以通过使用results()函数中的cooksCutoff参数来关闭这种过滤。

# Filter genes that have an extreme outlier
res_tableOE[which(is.na(res_tableOE$pvalue) & is.na(res_tableOE$padj) &res_tableOE$baseMean > 0),] %>% data.frame() %>% View()

如果一个基因包含一个有极端离群值的样本，那么p值和调整后的p值将被设置为NA。

平均归一化计数低的基因

DESeq2定义了一个较低的平均阈值，这是根据你的数据经验性确定的，在这个阈值中，可以通过减少被考虑进行多次检验的基因数量来增加显著基因的比例。这是基于这样一种观点，即计数非常低的基因不太可能看到显著差异，这通常是由于高离散。

在这里插入图片描述

在用户指定的值(alpha = 0.05)下，DESeq2根据平均计数过滤掉越来越大的部分基因，评估显著基因数量的变化，如上图所示。重要基因数量达到峰值的点是用于过滤经过多次检验的基因的低平均阈值。还有一个参数通过设置independentFiltering = F来关闭过滤。

# Filter genes below the low mean threshold
res_tableOE[which(!is.na(res_tableOE$pvalue) & is.na(res_tableOE$padj) & res_tableOE$baseMean > 0),] %>% data.frame() %>% View()

如果一个基因被独立过滤而过滤，那么只有调整后的p值被设置为NA。
注意:DESeq2将在默认情况下执行上述筛选;然而，其他DE工具，如EdgeR将不会。过滤是一个必要的步骤，即使您正在使用limma-voom和/或edgeR的拟似然方法（quasi-likelihood methods）。在使用其他工具时，一定要遵循预先过滤的步骤，如在Bioconductor上的用户指南中所述，因为它们通常表现得更好。

倍数变化（Fold change）

结果表中另一重要的列是log2FoldChange。由于有大量显著的基因列表，很难提取出有意义的生物学相关性。为了帮助增加强度，还可以添加倍数变化阈值。在设置这个值时要记住，我们使用的是log2 fold changes，因此log2FoldChange < 1的截断值将转换为实际的fold change为2。

另一种添加倍数变化阈值的方法:
results()函数有一个选项，可以使用lfcThrehsold参数添加倍数变化阈值。这种方法更具有统计学上的意义，当你想要一组基于特定倍数变化的更有信心的基因时，推荐使用这种方法。通过对大于指定绝对值的log2 倍数变化执行双尾检验，它实际上对期望的阈值执行统计检验。用户可以使用altHypothesis更改备择假设，并执行两个单尾检验。这是一种更保守的方法，所以期望检索到一组更小的基因!

结果表中报告的倍数变化计算公式为:

log2 (normalized_counts_group1 / normalized_counts_group2)

问题是，这些倍数变化估计并不完全准确，因为它们没有考虑到我们观察到的低读段计数的大离散度。为了解决这个问题，需要调整log2 fold changes。

更准确的LFC估计

为了产生更准确的log2倍数变化(LFC)估算，DESeq2考虑到当基因信息量较低时LFC估算趋近于零（向零收缩），这可能包括:

低计数
高离散值

LFC收缩使用来自所有基因的信息来产生更准确的估计。具体地说，所有基因的LFC估计值的分布(作为先验)被用来缩小信息少或高离散基因的LFC估计值，使其向更可能(更低)的LFC估计值靠拢。

在上图中，我们有一个使用两个基因的例子绿色基因和紫色基因。对于每个基因，我们绘制了两个不同小鼠品系(C57BL/6J和DBA/2J)中每个样本的表达值。两个基因在两个样本组中具有相同的均值，但绿色基因在组内变异较小，而紫色基因变异程度较高。对于组内变异较低的绿色基因，未缩小的LFC估计(绿色实线的顶点)与缩小的LFC估计(绿色虚线的顶点)非常相似。然而，由于高分散性，对紫色基因的LFC估计有很大不同。因此，即使两个基因可以有相似的标准化计数值，他们可以有不同程度的LFC收缩。注意，LFC的估计值向先验的方向缩小了(黑色实线)。
为了生成收缩的log2 fold change估计，你必须在使用lfcShrink()函数的results对象(我们将在下面创建)上运行额外的步骤。

## Save the unshrunken results to compare
res_tableOE_unshrunken <- res_tableOE# Apply fold change shrinkage
res_tableOE <- lfcShrink(dds, contrast=contrast_oe, type = "normal")

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

不同类型的收缩估计
根据DESeq2的不同版本，你使用的收缩估计默认方法会有所不同。可以通过在lfcShrink()函数中添加参数type来更改默认值，就像我们上面做的那样。指定“normal”方法可能会导致以下文本出现在您的控制台:

using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).Note that type='apeglm' and type='ashr' have shown to have less bias than type='normal'.
See ?lfcShrink for more details on shrinkage type, and the DESeq2 vignette.
Reference: https://doi.org/10.1093/bioinformatics/bty895

Normal是DESeq2中最初使用的方法。最近，作者实现了不同方法的收缩估计选项。已经证明，在大多数情况下，这些方法比我们默认使用的“正常”方法有更少的偏差。
但是，在使用这些方法时，您将需要指定coef，而不是使用contrast参数。使用contrasrt形成一个扩展的模型矩阵，平等地处理所有因素水平，并平均所有对因素水平之间的所有距离，以估计先验。使用coef，意味着只看模型矩阵的这一列(因此，这通常是相对于参考水平的一层)，并从这些系数的MLE分布中估计该系数的先验值。当使用coef时，收缩取决于选择哪一个水平作为参考。
例如，如果使用apeglm，你将希望运行以下代码:

res_table_apeglm <- lfcShrink(dds, coef="sampletype_MOV10_overexpression_vs_control", type="apeglm")

关于收缩的更多信息，DESeq2 vignette有一个关于收缩估算的扩展部分，这是非常有用的。

MA图

对探索我们结果有用的一个图是MA图。MA图显示了所有检测基因的归一化计数平均值与log2倍数变化。显着DE的基因被染色以便识别。这也是一个很好的方法来说明LFC收缩的影响。DESeq2包提供了一个简单的功能来生成一个MA plot。
让我们从未收缩（shrunken）的结果开始:

# MA plot using unshrunken fold changes
plotMA(res_tableOE_unshrunken, ylim=c(-2,2))

现在收缩后的结果是:

# MA plot using shrunken fold changes
plotMA(res_tableOE, ylim=c(-2,2))

在左边有未缩小的倍数变化值，可以看到大量的离散为低表达的基因。也就是说，许多低表达基因表现出非常高的倍数变化。收缩后，我们看到的倍数变化是估计值变小多了。
除了上面描述的比较，这个图允许我们评估倍数变化的幅度，以及它们相对于平均值的分布。一般来说，我们期望在所有的表达水平上看到显著的基因。

Excercise

MOV10 Differential Expression Analysis: Control versus Knockdown

Now that we have results for the overexpression results, do the same for the Control vs. Knockdown samples.

Create a contrast vector called contrast_kd.
Use contrast vector in the results() to extract a results table and store that to a variable called res_tableKD.
Shrink the LFC estimates using lfcShrink() and assign it back to res_tableKD.、

Answers

# 1. Create a contrast vector called contrast_kd.
contrast_kd <- c("sampletype", "MOV10_knockdown", "control")# 2. Use contrast vector in the results() to extract a results table and store that to a variable called res_tableKD.
res_tableKD <- results(dds, contrast=contrast_kd, alpha = 0.05)# 3. Shrink the LFC estimates using lfcShrink() and assign it back to res_tableKD.res_tableKD <- lfcShrink(dds, contrast=contrast_kd, type = "normal") #latest version of DESeq2
## OR
res_tableKD <- lfcShrink(dds, contrast=contrast_kd, res=res_tableKD) #older versions of DESeq2

哈佛大学——差异表达分析（十）Wald检验结果解读相关推荐

哈佛大学——差异表达分析（八）假设检验和多重检验校正
文章目录学习目标 DESeq2:模型拟合和假设检验广义线性模型(Generalized Linear Model) 假设检验 Wald 检验似然比检验(Likelihood ratio test ...
哈佛大学——差异表达分析（七）设计公式（Design formulas）
文章目录学习目标利用DESeq2进行差异表达分析运行DESeq2 设计公式(design formula) 复杂的设计 MOV10 差异表达分析学习目标使用DESeq2执行差异表达分析工作流 ...
哈佛大学——差异表达分析（九）DESeq2步骤描述
文章目录学习目标 DESeq2差异基因表达分析流程第一步:估计大小因子第二步:估计基因离散(gene-wise dispersion) 第三步:拟合曲线到基因的分散估计第四步:将基因离散估计值 ...
生物信息学入门使用 RNAseq counts数据进行差异表达分析（DEG）——edgeR 算法数据代码结果解读
差异表达分析通常作为根据基因表达矩阵进行生物信息学分析的第一步,有助于我们观察基因在不同样本中的表达差异,从而确定要研究的基因和表型之间的联系.常用的基因表达数据来自基因芯片或高通量测序.虽然矩阵看起 ...
MIT Technology Review 2020年“十大突破性技术”解读【中国科学基金】2020年第3期发布...
来源:国家自然科学基金委员会 MIT Technology Review 2020年"十大突破性技术"解读［编者按］ 2020年2月26日,MIT Technology Re ...
基因表达分析（上）- 差异表达分析
基因表达什么是基因表达,如下是来自于维基百科的解释: Gene expression is the process by which information from a gene is used ...
RNA-seq 详细教程： `DESeq2` 差异表达分析（7）
学习目标了解如何设计公式了解如何使用 DESeq2 执行差异表达分析 1. DE 分析差异表达分析工作流程的最后一步是将原始计数拟合到 NB 模型并对差异表达基因进行统计检验.在这一步中,我们本 ...
SCS【10】单细胞转录组之差异表达分析 (Monocle 3)
点击关注,桓峰基因桓峰基因公众号推出单细胞系列教程,有需要生信分析的老师可以联系我们!首选看下转录分析教程整理如下: Topic 6. 克隆进化之 Canopy Topic 7. 克隆进化之 Car ...
ballgown包进行基因差异表达分析
ballgown包可以读入Stringtie 的转录组组装及定量数据,进行基因差异表达分析. 1. 数据读入 # if (!requireNamespace("BiocManager&quo ...

哈佛大学——差异表达分析（十）Wald检验结果解读