DESeq2 — 基于负二项分布的 RNA-seq 差异表达分析包¶

# 在 R 中安装 DESeq2（通过 Bioconductor）
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# 安装 DESeq2（Bioconductor 3.22 对应 DESeq2 v1.50）
BiocManager::install("DESeq2")

# 安装常用辅助包
BiocManager::install(c("ggplot2", "pheatmap", "EnhancedVolcano"))
install.packages(c("ggplot2", "dplyr", "tibble"))

# 验证安装
library(DESeq2)
packageVersion("DESeq2")   # 查看版本号

library(DESeq2)
library(dplyr)

# 读取 count 矩阵（行=基因，列=样本）
# 数据来源：featureCounts/HTSeq/Salmon/Kallisto 的 count 输出
count_matrix <- read.csv("counts_matrix.csv",
                          row.names = 1,     # 第一列作行名（基因名）
                          check.names = FALSE) # 不转换列名格式

# 确保 count 矩阵为整数（DESeq2 要求）
count_matrix <- round(count_matrix)           # 四舍五入为整数

# 读取样本信息表（行=样本，列=分组信息等）
coldata <- data.frame(
    sample = colnames(count_matrix),          # 样本名
    condition = c("T2D", "T2D", "T2D",        # 条件分组（与 count 矩阵列顺序一致！）
                  "HC", "HC", "HC"),
    row.names = "sample"
)

# 关键检查：样本名顺序必须与 count 矩阵列顺序完全一致
all(rownames(coldata) == colnames(count_matrix))  # 必须返回 TRUE

# 构建 DESeqDataSet 对象
# countData：count 矩阵
# colData：样本信息表
# design：实验设计公式（~分组变量）
dds <- DESeqDataSetFromMatrix(
    countData = count_matrix,
    colData = coldata,
    design = ~ condition            # 比较 condition 分组
)

# 预过滤：去除几乎没有表达的基因（加速运算，非必须）
keep <- rowSums(counts(dds)) >= 10  # 保留总 count >= 10 的基因
dds <- dds[keep, ]

# 运行 DESeq2（核心函数，自动完成归一化+估算离散度+差异检验）
dds <- DESeq(dds)

# 查看比较组信息
resultsNames(dds)                   # 显示可用的对比关系

# 提取结果（默认 T2D vs HC）
# contrast：指定比较的分组
# alpha：FDR 显著性阈值
res <- results(dds,
               contrast = c("condition", "T2D", "HC"),  # 比较组设置
               alpha = 0.05)            # 显著性阈值

# 查看结果摘要
summary(res)                            # 显示上调、下调基因数量

# 按 adjusted p-value 排序
res_ordered <- res[order(res$padj), ]
head(res_ordered)                       # 查看最显著的基因

# 筛选差异显著基因（|log2FC| > 1 且 padj < 0.05）
deg <- subset(res_ordered,
              padj < 0.05 & abs(log2FoldChange) > 1)
nrow(deg)                               # 统计 DEG 数量

# 保存结果到文件
write.csv(as.data.frame(res_ordered),
          "DESeq2_results.csv",
          row.names = TRUE)

# 1. 火山图
library(EnhancedVolcano)
EnhancedVolcano(res,
    lab = rownames(res),               # 基因名标签
    x = "log2FoldChange",              # X 轴：log2FC
    y = "padj",                        # Y 轴：校正 p 值
    title = "T2D vs HC",
    pCutoff = 0.05,                    # p 值阈值线
    FCcutoff = 1.0)                    # FC 阈值线

# 2. MA 图（均值 vs 变化）
plotMA(res, ylim = c(-5, 5))           # 评估归一化效果

# 3. 热图（Top 50 DEG）
library(pheatmap)
vst_data <- vst(dds, blind = FALSE)   # 方差稳定化转换（用于可视化）
top_genes <- head(rownames(deg), 50)   # 取前 50 个 DEG
pheatmap(assay(vst_data)[top_genes, ],
         scale = "row",                # 行方向标准化
         cluster_cols = TRUE,          # 列（样本）聚类
         annotation_col = coldata)     # 标注分组信息

# 4. PCA 图（样本质量控制）
plotPCA(vst_data, intgroup = "condition") +
    ggplot2::theme_bw()

# 完整最简流程（6行核心代码）
library(DESeq2)
counts <- read.csv("counts.csv", row.names=1)
coldata <- data.frame(condition=c("T2D","T2D","HC","HC"), row.names=colnames(counts))
dds <- DESeqDataSetFromMatrix(counts, coldata, ~ condition)
dds <- DESeq(dds)
res <- results(dds, contrast=c("condition","T2D","HC"), alpha=0.05)
write.csv(as.data.frame(res), "results.csv")