引言

本系列[1]将开展全新的转录组分析专栏，主要针对使用DESeq2时可能出现的问题和方法进行展开。

提升速度与并行计算的思考

对于大多数分析任务，上述步骤的耗时通常不会超过30秒。然而，面对那些设计复杂且样本众多的实验（比如，涉及几十个系数和大约100个样本），用户可能需要比DESeq默认运行更快的计算速度。为此，提供两点建议：

1. 首先，通过设置参数 fitType="glmGamPoi"，用户可以利用由Constantin Ahlmann-Eltze开发的更高效的负二项广义线性模型（NB GLM）引擎。需要注意的是，在DESeq2中使用glmGamPoi时，必须指定测试类型为"LRT"，并定义一个简化的设计模型。
2. 其次，用户可以利用并行计算的优势。通过加载BiocParallel包，并设置参数parallel=TRUE和BPPARAM=MulticoreParam(4)，用户可以轻松实现DESeq、结果和lfcShrink的并行处理，例如，将任务分配到4个核心上。不过，在进行并行计算时，有几点建议需要考虑：首先，建议预先筛选掉那些在所有样本中计数都很低的基因，以避免无谓地将数据发送到子进程，因为这些基因的统计力低，最终也会被独立筛选掉；其次，增加更多核心往往会因为数据传输的开销而导致收益递减，因此建议从增加少量核心开始尝试。另外，对于resultsNames(dds)中列出的系数或对比结果的获取速度很快，无需进行并行处理。作为BPPARAM的另一种选择，用户可以在分析开始时注册核心，然后在调用函数时只需设置parallel=TRUE即可。

library("BiocParallel")
register(MulticoreParam(4))

P值与校正后的P值

可以根据最小的P值对结果表格进行排序。

resOrdered <- res[order(res$pvalue),]

可以使用汇总函数总结一些基本的计数。

summary(res)

image-20241123180220726

有多少个调整后的 p 值小于 0.1？

sum(res$padj < 0.1, na.rm=TRUE)

## [1] 1069

结果函数提供了多种参数，用于定制生成的结果表格。你可以通过查询 ?results 来获取这些参数的详细信息。需要注意的是，结果函数会自动根据每个基因的标准化计数均值进行独立筛选，以优化在给定的假发现率（FDR）阈值，即 α 下，拥有调整后 p 值低于该阈值的基因数量。独立筛选的更多细节将在后文讨论。默认情况下，参数 α 被设定为 0.1。如果调整后的 p 值阈值不是 0.1，那么应该将 α 设置为那个特定的值。

res05 <- results(dds, alpha=0.05)
summary(res05)

image-20241123180339717

sum(res05$padj < 0.05, na.rm=TRUE)

## [1] 853

独立假设权重分配

p值筛选概念的一个扩展是，通过对假设进行权重分配来提升检测能力。Bioconductor 提供了一个名为 IHW 的包，该包实现了独立假设权重（IHW）的方法。在此，展示了如何利用 IHW 对 DESeq2 的结果进行 p 值校正。

# (unevaluated code chunk)
library("IHW")
resIHW <- results(dds, filterFun=ihw)
summary(resIHW)
sum(resIHW$padj < 0.1, na.rm=TRUE)
metadata(resIHW)$ihwResult

高级用户须知，DESeq2 包计算出的所有数据均保存在 DESeqDataSet 或 DESeqResults 对象中，如何获取这些数据的讨论将在后文展开。

MA图

在 DESeq2 中，plotMA 函数用于展示在 DESeqDataSet 中所有样本的标准化计数均值基础上，由特定变量引起的 log2 倍数变化。如果调整后的 p 值小于 0.1，相应的点将以蓝色标出。超出视窗范围的点将以空心的向上或向下的三角形表示。

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

image-20241123180424090

可视化缩小的 log2 倍数变化的 MA 图更有用，它可以消除与低计数基因的 log2 倍数变化相关的噪声，而不需要任意的过滤阈值。

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

在执行了 plotMA 函数之后，用户可以利用 identify 函数通过点击图表来交互式地确定特定基因的行号。接着，用户可以通过保存这些索引值来检索基因的标识符。

idx <- identify(res$baseMean, res$log2FoldChange)

rownames(res)[idx]

替代的收缩估计方法

Love、Huber 和 Anders提出的调整对数倍数变化（log fold changes）基于一个以零为中心的正态先验分布，其尺度参数会根据数据进行适配。这种收缩后的对数倍数变化有助于进行排序和可视化，且无需对低计数基因设置任意的筛选阈值。然而，正态先验在某些数据集上可能会导致过度收缩。

• apeglm 是 apeglm 包中的自适应 t 分布先验收缩估计方法。从 1.28.0 版本起，它成为了默认的估计方法。
• ashr 是 ashr 包中的自适应收缩估计方法。DESeq2 利用 ashr 选项来适配一个正态分布的混合模型作为先验，设置 method="shrinkage"。
• normal 是 DESeq2 最初的收缩估计方法，它使用了一个自适应的正态分布作为先验。

在上述的 LFC 收缩代码示例中，明确指定了系数为 "
condition_treated_vs_untreated"。或者，也可以通过系数在 resultsNames(dds) 中的顺序来指定它，例如这里可以简单地使用 coef=2。

resultsNames(dds)

## [1] "Intercept"                      "condition_treated_vs_untreated"

# because we are interested in treated vs untreated, we set 'coef=2'
resNorm <- lfcShrink(dds, coef=2, type="normal")
resAsh <- lfcShrink(dds, coef=2, type="ashr")

par(mfrow=c(1,3), mar=c(4,4,2,1))
xlim <- c(1,1e5); ylim <- c(-3,3)
plotMA(resLFC, xlim=xlim, ylim=ylim, main="apeglm")
plotMA(resNorm, xlim=xlim, ylim=ylim, main="normal")
plotMA(resAsh, xlim=xlim, ylim=ylim, main="ashr")

image-20241123180701044

提示：已经优化了 apeglm 方法，现在它的运行时间与 normal 方法相近，例如，处理大约包含 10,000 个基因和 7 个样本的 pasilla 数据集大约需要 5 秒。如果需要快速估算 LFC 的收缩值，但又不需要后验标准差，可以将 apeMethod 设置为 "nbinomC"，这将大幅提升速度（约 10 倍），但会导致 lfcSE 列的值变为不可用（NA）。apeMethod 的另一个选项 "nbinomC*" 包含了随机启动，是一种快速方法的变体。

提示：如果数据中存在不希望的变异（例如，批次效应），建议进行校正。在 DESeq2 中，可以通过在设计中包含已知的批次变量，或者使用 sva 包中的 svaseq 函数/包或 RUVSeq 中的 RUV 函数来估计并校正这些不希望的变异。ashr 的开发者还提供了一种特别的方法，用于结合 ashr 来处理不希望的变异。

绘制计数图

检查单个基因在不同组中的读数计数有时也是有益的。plotCounts 函数可以简单地实现这一绘图功能，它根据估计的大小因子（或如果使用了这些，则是归一化因子）来归一化计数，并添加一个 1/2 的伪计数以支持对数刻度绘图。计数会根据 intgroup 中的变量进行分组，允许指定多个变量。在此，选择了上述结果表中 p 值最小的基因。你可以通过基因名称或数字索引来选择要绘制的基因。

plotCounts(dds, gene=which.min(res$padj), intgroup="condition")

对于自定义绘图，参数 returnData 指定该函数应仅返回用于使用 ggplot 绘图的 data.frame。

d <- plotCounts(dds, gene=which.min(res$padj), intgroup="condition", 
                returnData=TRUE)
library("ggplot2")
ggplot(d, aes(x=condition, y=count)) + 
  geom_point(position=position_jitter(w=0.1,h=0)) + 
  scale_y_log10(breaks=c(25,100,400))