deseq2-differential-expression

DESeq2 差异表达分析（R/Bioconductor）

概述

DESeq2 是 Bioconductor 的 R 包，用于分析 bulk RNA-seq 计数数据中的差异基因表达。它为每个基因拟合负二项广义线性模型，使用跨基因的经验贝叶斯收缩估计离散度参数，并使用 Wald 检验（两组比较）或似然比检验（复杂设计）检测差异表达。DESeq2 是 RNA-seq 差异表达分析领域的 R 金标准，并与 Bioconductor 原生集成，可无缝导入 Salmon（tximeta/tximport）、featureCounts 或 HTSeq 的数据。

适用场景

• 根据 bulk RNA-seq 计数数据，找出两种实验条件（处理组与对照组、疾病组与健康组）之间的差异表达基因
• 分析考虑批次效应或协变量的多因素设计（例如，~ batch + condition）
• 使用交互项（例如，时间 × 处理）或简化模型，通过似然比检验（LRT）检验复杂假设
• 通过 tximeta 或 tximport 导入 Salmon 假比对结果，以传播转录本层面的不确定性
• 使用 apeglm 进行 LFC 收缩，生成排序后的基因列表、火山图并开展下游通路分析
• 开展时间序列实验，或分析任何需要进行模型比较且水平数超过两个的设计
• 在 R/Bioconductor 生态系统中工作，需要与 SummarizedExperiment、clusterProfiler 或 EnhancedVolcano 集成
• 对于基于 Python 的流程，并需要使用相同的统计模型，请改用 pydeseq2-differential-expression
• 对于使用 TMM 归一化、准似然 F 检验或 TREAT 检验的负二项差异表达分析，请使用 edgeR
• 差异表达分析后，以及需要为 DESeq2 结果生成符合发表要求的图表时，请使用 omics-plotting SKILL
• 差异表达分析后，需要在通路层面解读结果时，请使用 gseapy-gene-enrichment SKILL

前置条件

• R 包：DESeq2（Bioconductor）、tximeta 或 tximport（导入 Salmon 数据）、apeglm（LFC 收缩）、pheatmap、ggplot2、EnhancedVolcano
• 数据要求：原始（未归一化）整数计数矩阵——行是基因、列是样本——以及列名与之匹配的样本元数据数据框。如果使用 Salmon：每个样本的 quant.sf 文件和转录本到基因的映射
• 环境：R ≥ 4.2、Bioconductor ≥ 3.16

# Install Bioconductor packages
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("DESeq2", "tximeta", "tximport", "apeglm",
                        "EnhancedVolcano"))
install.packages(c("pheatmap", "ggplot2", "dplyr"))

快速开始

使用计数矩阵完成两组比较，不到 20 行代码。

library(DESeq2)

# Load count matrix (genes x samples) and metadata
counts <- as.matrix(read.csv("counts.csv", row.names = 1))
coldata <- read.csv("metadata.csv", row.names = 1)
coldata$condition <- factor(coldata$condition)

# Build DESeqDataSet, run full pipeline
dds <- DESeqDataSetFromMatrix(countData = counts,
                               colData   = coldata,
                               design    = ~ condition)
dds <- dds[rowSums(counts(dds)) >= 10, ]   # pre-filter
dds <- DESeq(dds)

# Extract results (treated vs. control)
res <- results(dds, contrast = c("condition", "treated", "control"),
               alpha = 0.05)
summary(res)
# LFC shrinkage for visualization
res_shrunk <- lfcShrink(dds, contrast = c("condition", "treated", "control"),
                         type = "apeglm", coef = 2)

工作流

步骤 1：准备计数矩阵和样本元数据

根据基因 × 样本计数矩阵和 colData 数据框构建 DESeqDataSet。计数矩阵的列名必须与 colData 的行名一致。

library(DESeq2)
library(dplyr)

# Load raw count matrix (genes as rows, samples as columns)
counts <- as.matrix(read.csv("featureCounts_matrix.csv", row.names = 1))
# Or from featureCounts tab-delimited output — skip first comment line
# fc <- read.table("featurecounts_output.txt", header = TRUE, skip = 1, row.names = 1)
# counts <- as.matrix(fc[, 7:ncol(fc)])  # columns 7+ are sample counts

# Load sample metadata
coldata <- data.frame(
    condition = factor(c("control", "control", "control",
                          "treated", "treated", "treated")),
    batch     = factor(c("A", "A", "B", "A", "B", "B")),
    row.names = colnames(counts)
)

# Build DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = counts,
                               colData   = coldata,
                               design    = ~ condition)

cat("Samples:", ncol(dds), "\n")
cat("Genes:", nrow(dds), "\n")
cat("Condition levels:", levels(coldata$condition), "\n")

步骤 2：通过 tximeta 从 Salmon 导入数据

如果使用 Salmon 对 reads 进行定量，请使用 tximeta 导入转录本层面的估计值，并进行适当的偏移校正，以考虑转录本长度和 GC 偏差。

library(tximeta)

# Build a coldata with a "files" column pointing to quant.sf files
quant_dirs <- file.path("salmon_output", coldata$sample_id, "quant.sf")
coldata$files <- quant_dirs
coldata$names <- coldata$sample_id

# Import with tximeta — automatically fetches transcript metadata from Ensembl
se <- tximeta(coldata)           # SummarizedExperiment with transcript-level data

# Summarize to gene level (requires Ensembl or custom txdb)
gse <- summarizeToGene(se)       # gene-level SummarizedExperiment

# Build DESeqDataSet from the SummarizedExperiment
dds <- DESeqDataSet(gse, design = ~ condition)

cat("Gene-level SE dimensions:", dim(gse), "\n")
cat("Assay names:", assayNames(gse), "\n")
# Assay "counts" = estimated counts; "abundance" = TPM; "length" = eff. lengths

步骤 3：预过滤和质量控制

移除计数极低的基因，以提高统计功效并减少多重检验负担。使用方差稳定化后的计数进行 PCA，检查样本质量。

library(ggplot2)

# Pre-filter: keep genes with at least 10 reads total
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]
cat("Genes after pre-filtering:", nrow(dds), "\n")

# Variance-stabilizing transformation for QC visualization
vsd <- vst(dds, blind = TRUE)    # blind = TRUE for exploratory QC

# PCA plot
pca_data <- plotPCA(vsd, intgroup = c("condition", "batch"), returnData = TRUE)
percent_var <- round(100 * attr(pca_data, "percentVar"))

ggplot(pca_data, aes(PC1, PC2, color = condition, shape = batch)) +
    geom_point(size = 3) +
    xlab(paste0("PC1: ", percent_var[1], "% variance")) +
    ylab(paste0("PC2: ", percent_var[2], "% variance")) +
    ggtitle("PCA — Variance-Stabilized Counts") +
    theme_bw()

ggsave("pca_plot.pdf", width = 6, height = 5)
cat("Saved pca_plot.pdf\n")

步骤 4：运行 DESeq()——归一化、离散度估计和模型拟合

DESeq() 会依次执行三个步骤：（1）使用中位数比值法估计 size factor，（2）估计基因层面及收缩后的离散度，（3）拟合负二项 GLM 并计算 Wald 统计量。

# Run full DESeq2 pipeline
dds <- DESeq(dds)

# Inspect size factors (library normalization)
cat("Size factors:\n")
print(sizeFactors(dds))
# Expected range: 0.3–3.0; outliers indicate library quality issues

# Dispersion plot — should show decreasing dispersion trend with mean
plotDispEsts(dds, main = "Dispersion Estimates")

# Inspect fitted model coefficients
resultsNames(dds)
# [1] "Intercept" "condition_treated_vs_control"

步骤 5：提取结果并进行 FDR 校正

使用 results() 提取 Wald 检验统计量。为明确起见，请显式指定 contrast。独立过滤可最大限度地增加可检测基因的数量。

# Extract results for the contrast of interest
res <- results(dds,
               contrast          = c("condition", "treated", "control"),
               alpha             = 0.05,   # FDR threshold for independent filtering
               lfcThreshold      = 0,      # Test H0: |LFC| = 0 (Wald test)
               independentFilter = TRUE)   # Maximize power via adaptive filtering

# Summary: how many genes are up/down/NA
summary(res)
# out of N with nonzero total read count:
# adjusted p-value < 0.05: ...
# LFC > 0 (up): ...
# LFC < 0 (down): ...
# outliers [Cook's distance]: ...
# low counts [below independent filtering threshold]: ...

# Convert to data frame and sort by adjusted p-value
res_df <- as.data.frame(res)
res_df <- res_df[order(res_df$padj, na.last = TRUE), ]

cat("Significant genes (padj < 0.05):", sum(res_df$padj < 0.05, na.rm = TRUE), "\n")
cat("Upregulated:", sum(res_df$padj < 0.05 & res_df$log2FoldChange > 0, na.rm = TRUE), "\n")
cat("Downregulated:", sum(res_df$padj < 0.05 & res_df$log2FoldChange < 0, na.rm = TRUE), "\n")

write.csv(res_df, "deseq2_all_results.csv")

步骤 6：使用 apeglm 进行 LFC 收缩

使用自适应 t 先验（apeglm）收缩 log2 倍数变化。将收缩后的 LFC 用于 MA 图、火山图和基因排序——不要用于显著性判定（padj 来自未收缩的 Wald 检验）。

library(apeglm)

# apeglm shrinkage — specify coefficient name from resultsNames(dds)
res_shrunk <- lfcShrink(dds,
                         coef = "condition_treated_vs_control",
                         type = "apeglm",
                         res  = res)    # pass original results for speed

# Compare shrunk vs. unshrunk LFC
cat("Max |LFC| unshrunk:", max(abs(res$log2FoldChange), na.rm = TRUE), "\n")
cat("Max |LFC| shrunk:  ", max(abs(res_shrunk$log2FoldChange), na.rm = TRUE), "\n")

# Export shrunk results
res_shrunk_df <- as.data.frame(res_shrunk)
res_shrunk_df <- res_shrunk_df[order(res_shrunk_df$padj, na.last = TRUE), ]
write.csv(res_shrunk_df, "deseq2_shrunk_results.csv")

# Filter significant with effect size threshold
sig_genes <- res_shrunk_df[
    !is.na(res_shrunk_df$padj) &
    res_shrunk_df$padj < 0.05 &
    abs(res_shrunk_df$log2FoldChange) > 1, ]
cat("Significant (padj<0.05, |LFC|>1):", nrow(sig_genes), "\n")
write.csv(sig_genes, "deseq2_significant.csv")

步骤 7：可视化——MA 图、火山图和热图

计算与可视化分离（必需）： 不要在同一个代码块中同时进行 DESeq2 计算和绘图。步骤 6（或步骤 5）会保存结果 CSV；此步骤加载该 CSV 并仅生成图表。这样无需重新运行耗时的计算，即可重新生成图表或调整其样式。 绘图： 对于火山图，请在导出的 DEG CSV 上阅读 skills/data-visualization/omics-plotting/SKILL.md 并遵循其中的“Volcano”方案（→ figures/volcano_plot.pdf）。下方的热图仍使用原生 R（pheatmap）。

# Volcano: use the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "Volcano" recipe on the exported DEG CSV -> figures/volcano_plot.pdf

library(pheatmap)
# ---Heatmap of top 50 DE genes ---
top50 <- head(rownames(sig_genes[order(sig_genes$padj), ]), 50)
mat   <- assay(vsd)[top50, ]            # variance-stabilized expression
mat   <- mat - rowMeans(mat)            # center by gene mean

annotation_col <- data.frame(
    Condition = coldata$condition,
    Batch     = coldata$batch,
    row.names = colnames(mat)
)
# Use the code in the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) for publication-quality heatmap plots or use below code to generate a heatmap
pheatmap(mat,
         annotation_col  = annotation_col,
         cluster_rows    = TRUE,
         cluster_cols    = TRUE,
         show_rownames   = TRUE,
         fontsize_row    = 6,
         filename        = "figures/heatmap_top50.pdf",
         width           = 8,
         height          = 10)

步骤 8：多因素设计与交互项

适用于需要批次校正、配对样本或两个变量之间存在交互效应的设计。

# --- Batch-corrected two-group comparison ---
dds_batch <- DESeqDataSetFromMatrix(countData = counts,
                                     colData   = coldata,
                                     design    = ~ batch + condition)
dds_batch <- dds_batch[rowSums(counts(dds_batch)) >= 10, ]
dds_batch <- DESeq(dds_batch)

res_batch <- results(dds_batch,
                      contrast = c("condition", "treated", "control"),
                      alpha    = 0.05)
cat("Significant (batch-corrected):", sum(res_batch$padj < 0.05, na.rm = TRUE), "\n")

# --- Interaction design: test whether treatment effect differs between genotypes ---
# H0: the effect of treatment is the same in WT and KO
coldata$genotype <- factor(coldata$genotype)  # "WT" or "KO"

dds_int <- DESeqDataSetFromMatrix(
    countData = counts, colData = coldata,
    design = ~ genotype + condition + genotype:condition
)
dds_int <- DESeq(dds_int)
resultsNames(dds_int)
# Interaction coefficient: "genotypeKO.conditiontreated"

# Extract interaction term — genes where treatment effect differs by genotype
res_int <- results(dds_int, name = "genotypeKO.conditiontreated", alpha = 0.05)
cat("Genes with significant interaction:", sum(res_int$padj < 0.05, na.rm = TRUE), "\n")

# --- Likelihood Ratio Test for complex designs (e.g., time series) ---
# LRT compares full model vs. reduced model
dds_lrt <- DESeq(dds_int, test = "LRT", reduced = ~ genotype + condition)
res_lrt  <- results(dds_lrt)
cat("Genes with genotype:condition interaction (LRT):",
    sum(res_lrt$padj < 0.05, na.rm = TRUE), "\n")

关键参数

参数	函数	默认值	范围 / 选项	作用
design	DESeqDataSetFromMatrix()	（必需）	~ var 或 ~ cov + var	指定 GLM 公式；将批次/协变量放在感兴趣变量之前
contrast	results()	NULL	c(var, level1, level2) 或系数名称	定义比较；为清晰起见，始终显式指定
alpha	results()	0.1	0.01–0.10	用于独立过滤的 FDR 阈值；不会改变 padj 值
lfcThreshold	results()	0	0–2	检验 H0：
independentFilter	results()	TRUE	TRUE/FALSE	自适应平均计数过滤，以最大化可拒绝假设的数量
type	lfcShrink()	"apeglm"	"apeglm"、"ashr"、"normal"	收缩先验；apeglm 速度最快且校准效果最佳
coef	lfcShrink()	（必需）	resultsNames() 中的系数名称	要收缩的模型系数；必须与 resultsNames() 的输出一致
test	DESeq()	"Wald"	"Wald"、"LRT"	Wald：两组比较；LRT：复杂设计的模型比较
reduced	DESeq()	NULL	不包含待检验项的公式	LRT 必需；指定零模型
blind	vst() / rlog()	TRUE	TRUE/FALSE	QC/探索时设为 TRUE；使用离散度信息进行下游分析时设为 FALSE

关键概念

负二项模型与离散度

DESeq2 使用负二项分布（均值 = µ，离散度 = α）对每个基因的计数建模，其中均值 µ 通过对数线性模型取决于实验设计。其关键创新是离散度收缩：通过经验贝叶斯方法，将基因层面的离散度估计值收缩至所有基因的拟合趋势，从其他基因中借用信息，以稳定小样本量下的估计值。

尺寸因子标准化

DESeq2 使用比值中位数方法：对每个样本，计算每个基因的计数与该基因在所有样本中的几何平均值之比，然后取所有基因比值的中位数作为尺寸因子。此方法对差异表达基因具有稳健性，且不需要对差异表达基因所占比例作任何假设。

Wald 检验与似然比检验

• Wald 检验（默认）：检验单个系数 β = 0；速度快，适用于两组比较和多因素设计中的特定对比
• LRT：通过似然比将完整模型与 reduced 模型进行比较；适用于检验多水平因子（例如 5 个时间点）是否有任何水平产生影响，或将交互项作为一个整体进行检验

Cook 距离离群值检测

DESeq2 会计算每个基因在每个样本中的 Cook 距离，以标记潜在的表达离群值。如果任一样本的 Cook 距离超过阈值，结果中会自动将该基因替换为 NA。使用 refit = TRUE（默认值）时，会排除离群样本，并针对该基因重新拟合模型。

常用方案

配方：从一个拟合模型中提取多个对比

适用于比较多个处理组——只需拟合一次 DESeq()，即可提取多个对比。

# Fit once
dds_multi <- DESeqDataSetFromMatrix(counts, coldata, design = ~ condition)
dds_multi <- dds_multi[rowSums(counts(dds_multi)) >= 10, ]
dds_multi <- DESeq(dds_multi)

# Extract contrasts — all vs. "control"
contrasts_list <- list(
    A_vs_ctrl = c("condition", "treatment_A", "control"),
    B_vs_ctrl = c("condition", "treatment_B", "control"),
    C_vs_ctrl = c("condition", "treatment_C", "control")
)

results_list <- lapply(names(contrasts_list), function(nm) {
    res <- results(dds_multi, contrast = contrasts_list[[nm]], alpha = 0.05)
    res_shrunk <- lfcShrink(dds_multi, contrast = contrasts_list[[nm]],
                             type = "apeglm", res = res)
    write.csv(as.data.frame(res_shrunk), paste0("results_", nm, ".csv"))
    cat(nm, "— significant:", sum(res$padj < 0.05, na.rm = TRUE), "\n")
    as.data.frame(res_shrunk)
})
names(results_list) <- names(contrasts_list)

配方：配对样本设计

将受试者纳入阻断因子，以消除受试者内部的变异。

# Paired design: each subject has pre- and post-treatment samples
coldata$subject   <- factor(c("S1","S1","S2","S2","S3","S3"))
coldata$timepoint <- factor(c("pre","post","pre","post","pre","post"))

dds_paired <- DESeqDataSetFromMatrix(counts, coldata,
                                      design = ~ subject + timepoint)
dds_paired <- dds_paired[rowSums(counts(dds_paired)) >= 10, ]
dds_paired <- DESeq(dds_paired)

res_paired <- results(dds_paired,
                       contrast = c("timepoint", "post", "pre"),
                       alpha    = 0.05)
cat("Significant (paired):", sum(res_paired$padj < 0.05, na.rm = TRUE), "\n")

配方：使用收缩后的 LFC × -log10(pvalue) 为 GSEA 基因排序

创建预排序基因列表，供 GSEA 或 gseapy 的 prerank 函数使用。

# Build ranking metric: signed -log10(pvalue) with LFC direction
res_rank <- as.data.frame(res_shrunk)
res_rank <- res_rank[!is.na(res_rank$pvalue) & !is.na(res_rank$log2FoldChange), ]
res_rank$rank_metric <- sign(res_rank$log2FoldChange) *
                         (-log10(res_rank$pvalue + 1e-300))
res_rank <- res_rank[order(res_rank$rank_metric, decreasing = TRUE), ]

# Export as two-column TSv for fgsea or gseapy prerank
write.table(data.frame(gene  = rownames(res_rank),
                        score = res_rank$rank_metric),
            "gsea_ranked_list.rnk",
            sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE)
cat("Ranked gene list saved:", nrow(res_rank), "genes\n")

配方：保存并重新加载已拟合的 DESeq2 对象

避免在后续会话中重新运行耗时的 DESeq() 步骤。

# Save fitted DDS object
saveRDS(dds, "dds_fitted.rds")
cat("Saved fitted DESeqDataSet to dds_fitted.rds\n")

# Reload in a new session
dds_loaded <- readRDS("dds_fitted.rds")
# Continue with results(), lfcShrink(), etc.
res_reloaded <- results(dds_loaded, contrast = c("condition", "treated", "control"))

预期输出

文件	描述
deseq2_all_results.csv	完整结果表：包含所有基因的 baseMean、log2FoldChange、lfcSE、stat、pvalue、padj
deseq2_shrunk_results.csv	使用 apeglm 收缩后的 LFC 结果；用于排序和可视化
deseq2_significant.csv	筛选后的结果：padj < 0.05 且
pca_plot.pdf	方差稳定化计数的 PCA；用于检查批次结构和离群值
dispersion_plot.pdf	基因级离散度与拟合离散度；应显示紧密分布在趋势线附近的点云
figures/volcano_plot.pdf	标注了显著性和 LFC 阈值的火山图
figures/heatmap_top50.pdf	前 50 个差异表达基因的层次聚类热图
gsea_ranked_list.rnk	用于通路富集分析（fgsea/gseapy）的预排序基因列表
dds_fitted.rds	序列化的 DESeqDataSet，用于检查点保存和恢复

故障排除

问题	原因	解决方案
Error: design contains one or more variables with all samples having the same value	设计因子只有一个水平	检查 table(coldata$condition)；确保存在所有因子水平
Model matrix not full rank	设计存在混杂（例如，批次与条件完全相关）	运行 table(coldata$batch, coldata$condition)；移除混杂变量或收集更多样本
所有或大多数 padj = NA	基因因独立过滤（平均计数较低）或 Cook's 离群值而被标记	检查 summary(res) 中被过滤基因的数量；放宽 alpha 或核查数据质量
生物学效应明显，但没有显著基因	实验效能不足或组内变异较大	确认每组 n ≥ 3；检查 PCA 中是否存在离群值；查看 p 值直方图（应在 0 附近出现峰值）
Error in lfcShrink: 未找到系数	coef 名称与 resultsNames(dds) 不匹配	运行 resultsNames(dds)，并将准确的系数字符串复制到 coef = 中
apeglm 因收敛警告而失败	LFC 估计值极端（数据稀疏或完全分离）	切换为 type = "ashr"，它对这类情况更稳健
size factor 非常大（> 5）	文库大小差异极大，或误将归一化计数作为输入	检查 colSums(counts(dds))；确保输入是来自比对器/计数器的原始整数计数
离散度图显示远离趋势线的离群点云	文库质量差或失败；元数据分组不正确	检查 PCA；移除离群样本；确认 design 公式与生物学分组相符

参考资料

• DESeq2 Bioconductor 页面 — 软件包主页、教程和参考手册
• DESeq2 教程 — 全面的示例，涵盖所有设计类型
• Love 等（2014），Genome Biology — DESeq2 原始论文：RNA-seq 数据中倍数变化和离散度的稳健估计
• Zhu 等（2019），Bioinformatics — 用于 DESeq2 的 apeglm LFC 收缩
• tximeta Bioconductor 页面 — 导入带有转录本级元数据的 Salmon 结果
• EnhancedVolcano Bioconductor 页面 — 用于差异表达结果的出版级火山图