edgeR — 基于负二项模型的数字基因表达差异分析包¶

一句话说明¶

edgeR 是 Bioconductor 经典 RNA-seq 差异表达包，采用精确统计检验和广义线性模型，特别擅长处理生物学重复少（甚至无重复）的实验设计（当前版本 v4+）。

安装与配置¶

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

BiocManager::install("edgeR")            # 当前稳定版 v4+

# 安装辅助包
install.packages(c("ggplot2", "limma"))
BiocManager::install("limma")

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

核心用法¶

准备数据并构建 DGEList 对象¶

library(edgeR)

# 读取 count 矩阵（行=基因，列=样本）
counts <- read.table("counts_matrix.txt",
                     header = TRUE,        # 第一行为列名
                     row.names = 1,        # 第一列为基因名
                     sep = "\t")           # 制表符分隔

# 定义分组信息
group <- factor(c("T2D", "T2D", "T2D", "HC", "HC", "HC"))

# 构建 DGEList 对象（edgeR 的核心数据结构）
dge <- DGEList(
    counts = counts,                       # count 矩阵
    group = group                          # 样本分组
)

# 过滤低表达基因（推荐：每组至少有一半样本 CPM > 1）
keep <- filterByExpr(dge, group = group)   # 自动计算过滤阈值
dge <- dge[keep, , keep.lib.sizes = FALSE] # 保留满足条件的基因
cat("过滤后基因数:", nrow(dge$counts), "\n")

TMM 归一化¶

# TMM（Trimmed Mean of M-values）归一化
# 是 edgeR 推荐的归一化方法，校正测序深度和 RNA 组成偏差
dge <- calcNormFactors(dge,
                       method = "TMM")     # 归一化方法（TMM/TMMwsp/RLE/upperquartile）

# 查看归一化因子（正常值应接近 1.0）
dge$samples$norm.factors                   # 显示每个样本的归一化因子

方法一：精确检验（适合单因素简单比较）¶

# 估算离散度
# 单个离散度值（common dispersion）
dge <- estimateDisp(dge)

# 精确检验（类似 Fisher 精确检验）
et <- exactTest(dge, pair = c("HC", "T2D"))  # 对比：T2D vs HC
topTags(et, n = 20)                           # 查看前 20 个显著基因

# 提取完整结果
result_et <- topTags(et, n = Inf)$table       # n=Inf 提取所有基因
result_et_sig <- subset(result_et,
                         FDR < 0.05 & abs(logFC) > 1) # 筛选 DEG

方法二：GLM 拟似然（推荐，适合复杂实验设计）¶

# 设计矩阵（支持多因素、协变量、批次效应校正）
design <- model.matrix(~ 0 + group)           # ~ 0 表示不含截距（便于对比设置）
colnames(design) <- levels(group)              # 列名设为分组名

# 估算离散度（GLM 模式）
dge <- estimateDisp(dge, design)               # 同时估算 common + trended + tagwise

# 拟似然 F 检验（QL F-test，比 LRT 更保守稳健）
fit <- glmQLFit(dge, design)                   # 拟合 GLM 模型

# 设置对比（T2D vs HC）
contrast <- makeContrasts(T2D - HC,
                           levels = design)    # 定义对比矩阵

# 执行差异检验
qlf <- glmQLFTest(fit, contrast = contrast)
result_qlf <- topTags(qlf, n = Inf)$table      # 提取全部结果

# 筛选显著 DEG
deg_qlf <- subset(result_qlf,
                   FDR < 0.05 & abs(logFC) > 1)
cat("显著差异基因数:", nrow(deg_qlf), "\n")

# 保存结果
write.csv(result_qlf, "edgeR_results.csv")

参数详解¶

函数/参数	说明	默认值
`filterByExpr`	自动过滤低表达基因	推荐使用
`calcNormFactors`	计算 TMM 归一化因子	method="TMM"
`estimateDisp`	估算基因离散度	必须
`glmQLFit`	拟似然 F 检验拟合	推荐
`glmQLFTest`	差异检验（QL）	推荐
`exactTest`	精确检验（单因素）	简单情况
`topTags`	提取排序结果	n=10

可视化¶

# 1. MD 图（Mean-Difference，类似 MA 图）
plotMD(qlf, main = "T2D vs HC")
abline(h = c(-1, 1), col = "blue", lty = 2) # 添加 FC 阈值线

# 2. BCV 图（生物变异系数）
plotBCV(dge)                                  # 展示离散度估算情况

# 3. MDS 图（样本聚类，类似 PCA）
plotMDS(dge,
        col = ifelse(group == "T2D", "red", "blue"),  # 按组着色
        main = "MDS Plot")
legend("topright", levels(group), col = c("red", "blue"), pch = 1)

# 4. 火山图
plot(result_qlf$logFC,
     -log10(result_qlf$FDR),
     xlab = "log2 Fold Change",
     ylab = "-log10(FDR)",
     pch = 20, cex = 0.5,
     col = ifelse(result_qlf$FDR < 0.05 & abs(result_qlf$logFC) > 1,
                  "red", "grey"))

常见报错与解决¶

报错信息	原因	解决方法
`No replicates`	无生物学重复	使用 `estimateDisp(method="common")`
`all zero counts`	存在全零基因	用 `filterByExpr` 过滤
`NaN in common dispersion`	数据问题	检查数据是否有负数或非整数
`glmFit failed`	设计矩阵秩亏	检查 design 矩阵，去除线性相关列
结果 logFC 过大	某样本中有零 count	添加 `prior.count=2` 平滑处理

速查表¶

# 精确检验（最简流程）
library(edgeR)
dge <- DGEList(counts=counts, group=group)
dge <- dge[filterByExpr(dge),]
dge <- calcNormFactors(dge)
dge <- estimateDisp(dge)
et <- exactTest(dge, pair=c("HC","T2D"))
res <- topTags(et, n=Inf)$table
write.csv(res, "edgeR_results.csv")

# GLM 拟似然检验（推荐）
design <- model.matrix(~0+group)
dge <- calcNormFactors(estimateDisp(dge, design))
fit <- glmQLFit(dge, design)
res <- topTags(glmQLFTest(fit, contrast=makeContrasts(T2D-HC, levels=design)), n=Inf)