DESeq2 差异表达分析怎么做对

跟完这篇，您能用一个原始 counts 矩阵和一张样本信息表，跑通 DESeq2 的标准三步，用两张诊断图确认结果可信，读懂结果表的每一列，定出一套站得住脚的差异基因筛选阈值，并绕开三类最常见的报错。

先备齐两样东西

DESeq2 的输入只有两样，缺一不可，且必须严丝合缝地对上：

表达矩阵 counts：行是基因、列是样本，每个格子是原始 read counts（非负整数）。
样本信息表 colData：行是样本、列是分组等变量（如 condition、batch）。它的行名必须和 counts 的列名一致。

library(DESeq2)

dim(counts)          # 基因数 × 样本数
head(colData)        # 每个样本属于哪一组
all(colnames(counts) == rownames(colData))   # 必须为 TRUE

最后一行返回 TRUE，才说明样本对齐了。这一步偷懒，后面所有结果都可能张冠李戴。矩阵怎么整理成这个形状，见 /tutorials/matrix-prep。

为什么必须用原始 counts

这是最容易踩、也最伤结论的一步。DESeq2 只接受原始 counts，不能喂 TPM、FPKM、RPKM，也不能喂已经归一化或取过对数的值。

原因在于它的统计模型：DESeq2 假设每个基因的 counts 服从负二项分布，并在内部自己估计文库大小因子（size factor）来做归一化。这个模型依赖「计数」这一事实——均值和方差之间有特定关系（计数越大、方差越大）。一旦您先把数据变成了 TPM 或对数值，这层均值—方差关系就被破坏了，内部归一化也会被重复施加，估计出来的差异和 p 值都不再可信。

上游工具	得到的值	怎么进 DESeq2
featureCounts、HTSeq	整数 counts	直接用 `DESeqDataSetFromMatrix()`
Salmon、kallisto	估计 counts（可能非整数）	用 `tximport` 导入，再走 `DESeqDataSetFromTximport()`

Salmon、kallisto 的估计 counts 不要自己四舍五入——tximport 为这种输入准备了专门的通路，还能顺带把转录本汇总到基因层面。各种表达量单位到底差在哪、该喂哪种，见 /docs/formats。

三步：建 dds → DESeq() → results()

标准流程就三步，记住顺序即可。

第一步，建对象并写清设计公式。 设计公式 ~ condition 告诉 DESeq2 您想比较的分组变量。

dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData   = colData,
  design    = ~ condition
)

# 明确指定对照组(参照水平),否则按字母顺序默认——这会决定倍数变化的方向
dds$condition <- relevel(dds$condition, ref = "control")

# 可选:预过滤极低表达基因,加快运算、减少多重检验负担
keep <- rowSums(counts(dds)) >= 10
dds  <- dds[keep, ]

第二步，一行跑完整个流程。 DESeq() 把估计 size factor、估计离散度、拟合负二项模型、做 Wald 检验这几件事一次做完。

dds <- DESeq(dds)

第三步，取出结果。 用 contrast 明确写清「谁比谁」，方向不会含糊。

res <- results(dds, contrast = c("condition", "treated", "control"))
summary(res)

# 常按校正后 p 值排序,方便看最显著的基因
res <- res[order(res$padj), ]
head(res)

补充一点：log2FoldChange 的原始估计在低表达基因上噪声很大。若要画火山图或做基因排序，惯例是再做一次收缩（shrinkage），让不可靠的大倍数变化回归：

resultsNames(dds)    # 先查可用的 coef 名字,下面这行要照抄

# type 可选 "apeglm" / "ashr" / "normal",apeglm 需另装同名包
resLFC <- lfcShrink(dds, coef = "condition_treated_vs_control", type = "apeglm")

收缩只影响 log2FoldChange 的展示，不改变 p 值与显著性判断。

跑完先看两张诊断图

在读结果表之前，先花两分钟确认这次分析没出岔子。两张图基本够用。

一、PCA 看样本分不分得开。 先做方差稳定变换，再看前两个主成分：

# vst 做方差稳定变换,只为可视化,不参与差异检验
vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = "condition")

怎么读这张图：

同组样本聚在一起、两组明显分开——理想情况，可以放心往下走。
有个别样本孤零零地远离同组——可能是离群样本，回去查建库或测序质量。
样本不按 condition 聚、反而按 batch 聚——说明有批次效应，设计公式里要加上 batch（见下文「常见报错」）。

二、MA 图看整体形态。 横轴是平均表达、纵轴是 log2FoldChange，padj 显著的点会高亮：

plotMA(resLFC, ylim = c(-3, 3))   # 用收缩后的结果,低表达处的散点会收敛

正常的 MA 图：点云大致对称地围绕 0 线，低表达（左侧）散得开、高表达（右侧）收得紧。如果整片点云系统性地偏向一侧，往往提示归一化或样本组成有问题，值得回头查。

结果表怎么读

results() 返回一张表，每行一个基因。逐列看懂它，比记住任何阈值都重要。

列	含义	怎么用
baseMean	所有样本归一化后表达量的平均	反映该基因整体表达高低；太低的基因倍数变化通常不稳
log2FoldChange	log2 尺度的倍数变化，正值上调、负值下调	log2FC=1 即 2 倍，=−1 即降到一半；方向由 contrast 决定
lfcSE	log2FoldChange 的标准误	越小说明这个倍数估计越稳
stat	Wald 检验统计量（≈ log2FC / lfcSE）	一般不直接看，是 pvalue 的来源
pvalue	原始 p 值	不要直接拿它当阈值——没做多重检验校正
padj	BH 法校正后的 p 值（即 FDR）	筛差异基因就看这一列

两个容易困惑的点：

padj 出现 NA 是正常的。 DESeq2 会对表达过低、几乎不可能显著的基因做独立过滤（independent filtering），直接把它们的 padj 记为 NA，从而减少多重检验的分母、提高检出力。另有个别基因因存在极端离群值，pvalue 也会是 NA。筛选时把 NA 当作「未通过」处理即可。
看 padj、不看 pvalue。 一次差异分析动辄检验上万个基因，不校正必然混入大量假阳性。

阈值怎么定

最常被问：显著性和倍数变化的线划在哪？一套广泛使用的默认是：

# 先把 padj 的 NA 排掉,再按两个条件筛
res_clean <- res[!is.na(res$padj), ]
sig <- subset(res_clean, padj < 0.05 & abs(log2FoldChange) > 1)
nrow(sig)   # 筛出多少个差异基因

即 padj < 0.05（FDR 控制在 5%）且 |log2FC| > 1（变化至少 2 倍）。

但请把它当常用惯例，而非铁律——它是起点，要按您的数据和问题调整：

padj 的线：0.05 最常见，也有人按数据的信噪比收紧到 0.01 或放宽到 0.1。关键是您筛的是 FDR，报告时要写清用了哪条线。
|log2FC| 的线：1（2 倍）是习惯值，并无生物学上的绝对依据。差异普遍偏小的体系（如某些临床样本）可放宽到 0.585（1.5 倍）；差异很大时也可收紧。
别只靠 pvalue 或 raw p<0.05 筛：那不是显著性阈值，会放进大量假阳性。
倍数变化不能单独用：一个 log2FC 很大但 padj 不显著的基因，往往只是低表达噪声，不能算差异基因。

拿到 sig 之后，常见的下一步是画火山图看整体分布（见 /tutorials/volcano），或把差异基因送去做富集分析（见 /tutorials/enrichment）。如果您对该用哪套阈值、甚至该不该用倍数变化门槛都拿不准，可以先看 /tutorials/choose-stats 里对不同场景的取舍。

常见报错

三类报错占了绝大多数，对症即可。

1. counts 不是整数。 报错提示大意是 assay 中存在非整数值。多半是您喂了 TPM/FPKM，或用了 Salmon/kallisto 的估计 counts。

若本就是原始整数 counts 却报错，检查是不是被读成了小数/字符型：counts <- round(as.matrix(counts))（仅在确认是真整数被误存时用）。
若来自 Salmon/kallisto，不要手动取整，改用 tximport + DESeqDataSetFromTximport()。

2. 样本名对不上。 counts 的列名与 colData 的行名不一致、或顺序不同，建对象时会报错或悄悄错位。

# 让 colData 的行严格按 counts 的列排序
colData <- colData[colnames(counts), ]
all(colnames(counts) == rownames(colData))   # 应为 TRUE

3. 设计公式写错。 几种典型：

分组变量被当成了数值。用 colData$condition <- factor(colData$condition) 显式转成因子。
忘了设参照水平，倍数变化方向和预期相反。用 relevel() 指定 ref，或在 results() 里用 contrast 写清。
想校正批次时把变量顺序写反了。关注变量放最后：design = ~ batch + condition，这样 results() 默认给的是校正批次后的 condition 效应。
报「model matrix is not full rank（模型矩阵不满秩）」：通常是两个变量完全共线（比如每个 batch 恰好只对应一种 condition），此时无法区分二者的效应，需要重新设计分组或去掉冗余变量。

在平台上跑一遍

不想配 R 环境、也不想跟报错缠斗，可以在百沐一下直接做这一步：上传原始 counts 矩阵和样本分组表，用自然语言说清「哪组比哪组、对照是谁、阈值怎么定」，平台会零代码替您跑完建 dds、DESeq()、results() 这三步，并把诊断图、结果表和火山图一并给出。您要做的只是核对分组对不对、阈值合不合理。

从上传到出图的完整流程见 /docs/omics-agent，或直接开始：https://app.baimuyixia.com