RNA-seq（3）：用 DESeq2 做差异表达分析——以 airway 数据为例

以解牛之法析生信，观微雀之形览科研。

今天我们一起来学习使用DESeq2进行转录组的差异分析。

三、DEseq2差异分析及对比

使用数据：airway
工具：DESeq2
相关R包：DESeq2（分析）、ggplot2（绘图）

1. 工具介绍：DESeq2

DESeq2 是 RNA-seq 差异表达分析 最主流的 R 包之一（发表于 2014 年，已被引用数万次），专门用于基于负二项分布模型的计数数据差异分析。它由 Bioconductor 项目维护，是转录组分析的标准工具。

输入数据：原始整数计数矩阵（不能是 TPM/FPKM），行=基因，列=样本

标准化方法：median-of-ratios 方法，通过计算每个基因的相对计数因子来校正文库大小差异，对异常值更稳健

差异检验：基于 Wald 检验 或 似然比检验（LRT），自动估计基因离散度

处理小样本：使用 经验贝叶斯方法 将基因间离散度信息“收缩”（shrinkage），使小样本（如 n=3）也能获得可靠结果

log2FC 收缩：提供 apeglm 和 ashr 两种收缩方法，对低计数基因的 log2FC 进行适度压缩，提高结果稳定性

DESeq2说明文档

Bioconductor 手册：https://bioconductor.org/packages/devel/bioc/manuals/DESeq2/man/DESeq2.pdf

github仓库：https://github.com/thelovelab/DESeq2

使用说明：

# 如何使用
DESeq(
object,
test = c("Wald", "LRT"),
fitType = c("parametric", "local", "mean", "glmGamPoi"),
sfType = c("ratio", "poscounts", "iterate"),
betaPrior,
full = design(object),
reduced,
quiet = FALSE,
minReplicatesForReplace = 7,
modelMatrixType,
useT = FALSE,
minmu = if (fitType == "glmGamPoi") 1e-06 else 0.5,
parallel = FALSE,
BPPARAM = bpparam()
)

2. 下载及加载包

# # 第一次安装可以选中ctrl+shift+c取消注释，运行安装
# 安装biocManager
# if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
# # 差异分析核心DESeq2包
# BiocManager::install(c("DESeq2", "apeglm"), update = FALSE)
# # 数据处理
# BiocManager::install(c("tidyverse", "tibble"), update = FALSE)
# # 绘图包
# BiocManager::install(c("ggplot2", "EnhancedVolcano", "ggrepel", "pheatmap", "ggpubr", "RColorBrewer"), update = FALSE)

# 1.2.3 加载包，确定已经安装了包后再加载
# 如果报错或者提示没有安装，返回上一步
library(DESeq2)
library(apeglm)
library(tidyverse)
library(ggplot2)
library(EnhancedVolcano)
library(ggrepel)
library(pheatmap)
library(ggpubr)

3. 示例数据加载：

airway数据：

if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
options(BioC_mirror = "https://mirrors.ustc.edu.cn/bioc/") # 切换中科大镜像下载
BiocManager::install("airway") #下载airway数据，2014年发表

library(airway) #加载包
data(airway) #加载airway数据集
airway #查看数据集是否加载成功

Bioconductor 页面：https://bioconductor.org/packages/release/data/experiment/html/airway.html

原始 GEO 项目：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778

airway 数据集来源于 GEO 项目 GSE52778，是该项目的子集（包含 4 个细胞系各 1 对处理/对照样本，共 8 个样本），数据处理方式与原始研究一致。

注：如果加载不成功，也可以手动下载，但是此步需要进一步处理才可使用。
地址：https://bioconductor.org/packages/release/data/experiment/html/airway.html
或者使用下载的GSE52778，但是也需要适当的处理才可使用。

4. 标准的 DESeq2 差异分析流程

可从airway中：获取原始整数技术矩阵并构建分组信息——过滤——构建 DESeq2 对象——运行 DESeq2，进行差异基因的比对。

# 标准的分析流程
# 2.1 提取数据
library(airway)
data(airway)

count <- assay(airway) # 原始整数计数矩阵
col_data <- colData(airway) # 样本注释（含 cell 和 dex）

# 构建分组信息（也可以直接从 col_data 提取）
group <- data.frame(
sample = colnames(count),
group = col_data$dex, # "trt" / "untrt"
cell = col_data$cell # 细胞系信息，用于校正批次
)

# 确认样本顺序一致（必须为 TRUE）
stopifnot(all(colnames(count) == group$sample))

# 2.2 过滤低表达基因（RNA-seq 标准操作，提高统计功效）
# 保留至少在 3 个样本中 count >= 10 的基因（可根据实际情况调整阈值）
keep <- rowSums(count >= 10) >= 3
count_filter <- count[keep, ]
cat("过滤前基因数：", nrow(count), "\\n")
cat("过滤后基因数：", nrow(count_filter), "\\n")

# 2.3 构建 DESeq2 对象（配对设计，校正细胞系效应）
dds <- DESeqDataSetFromMatrix(
countData = count_filter, # 使用过滤后的矩阵
colData = group,
design = ~ cell + group # 加入 cell 作为协变量
)

# 明确指定对照组（确保 log2FC 方向正确）
dds$group <- relevel(dds$group, ref = "untrt")

# 2.4 运行 DESeq2
dds <- DESeq(dds)
# 查看可用的系数名（用于后续对比）
resultsNames(dds)

可以看到过滤前后基因数为，过滤前基因数： 63677 ，过滤后基因数： 16596 。此步完成了差异分析，核心包为DESeq2，下步进行差异分析的提取。

5. 提取差异基因

提取差异基因，并且分别保存全部的分析结果和显著的差异基因。

# 2.5 提取两组比较的结果
# contrast 参数：c(分组列名, 实验组, 对照组)
res <- results(dds, contrast = c("group", "trt", "untrt"))

# 按 p 值排序（最显著的在前）
res <- res[order(res$pvalue), ]

# 转为数据框，并添加基因名列
res_df <- as.data.frame(res)
res_df$gene <- rownames(res_df)

# 保存全部基因的差异结果
write.csv(res_df, file = "DESeq2_all_gene_result.csv", row.names = FALSE)

# 2.6 筛选显著差异基因（阈值：|log2FC| > 1 且 padj < 0.05）
res_df$change <- "Stable"
res_df$change[res_df$log2FoldChange > 1 & res_df$padj < 0.05] <- "Up"
res_df$change[res_df$log2FoldChange < –1 & res_df$padj < 0.05] <- "Down"

# 统计上下调基因数量
table(res_df$change)
cat("显著上调基因数：", sum(res_df$change == "Up"), "\\n")
cat("显著下调基因数：", sum(res_df$change == "Down"), "\\n")

# 2.7 （可选）保存显著差异基因子集
sig_genes <- res_df[res_df$change != "Stable", ]
write.csv(sig_genes, file = "DESeq2_sig_genes.csv", row.names = FALSE)

# 打印前几个显著基因
head(res_df[res_df$change != "Stable", ])

可以看到 有显著上调基因511，下调基因489。

运行到这里我们完成了：提取差异基因，并且保存全部的分析结果和显著的差异基因结果。

后续绘制火山图、富集分析，就可直接使用这数据啦！

今天的生信知识分享就到这里啦！咱们下期见！

以解牛之法析生信，观微雀之形览科研。