limma — 线性模型微阵列/RNA-seq 差异分析包¶

一句话说明¶

limma 用线性模型和经验贝叶斯方差估算进行差异表达分析，是微阵列数据的黄金标准，结合 voom 转换后同样适用于 RNA-seq count 数据，并支持复杂多因素实验设计。

安装与配置¶

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

BiocManager::install("limma")            # 安装 limma

# 安装辅助包
install.packages(c("ggplot2", "statmod"))

# 验证安装
library(limma)
packageVersion("limma")                  # 查看当前版本

核心用法¶

limma-voom 流程（RNA-seq count 数据推荐）¶

library(limma)
library(edgeR)   # 需要 edgeR 做预处理（filterByExpr 和归一化）

# 1. 准备数据
counts <- read.table("counts.txt", header=TRUE, row.names=1, sep="\t")
group <- factor(c("T2D","T2D","T2D","HC","HC","HC"))

# 2. 过滤低表达基因（同 edgeR 流程）
dge <- DGEList(counts=counts, group=group)
keep <- filterByExpr(dge, group=group)
dge <- dge[keep,, keep.lib.sizes=FALSE]
dge <- calcNormFactors(dge, method="TMM")  # TMM 归一化

# 3. 构建设计矩阵（实验设计）
design <- model.matrix(~ 0 + group)       # 无截距模式
colnames(design) <- levels(group)          # 列名设为分组名
print(design)                              # 确认设计矩阵正确

# 4. voom 转换（核心：将 count 数据转为适合线性模型的连续值）
# voom 通过估算均值-方差关系，为每个观测值赋予精度权重
v <- voom(dge,
          design,                          # 实验设计矩阵
          plot = TRUE)                     # 绘制 mean-variance 趋势图（质控用）

线性模型拟合与差异检验¶

# 5. 拟合线性模型
fit <- lmFit(v, design)                    # 最小二乘拟合

# 6. 设置对比矩阵（定义比较关系）
contrast.matrix <- makeContrasts(
    T2D_vs_HC = T2D - HC,                  # 对比1：T2D vs HC
    levels = design                         # 参照设计矩阵
)

# 可以定义多个对比：
# contrast.matrix <- makeContrasts(
#     T2D_vs_HC = T2D - HC,
#     T1D_vs_HC = T1D - HC,
#     levels = design)

# 7. 应用对比并计算
fit2 <- contrasts.fit(fit, contrast.matrix) # 对比拟合

# 8. 经验贝叶斯方差压缩（核心：借用多基因信息提高稳健性）
fit2 <- eBayes(fit2)

# 9. 提取差异结果
result <- topTable(fit2,
                   coef = "T2D_vs_HC",     # 指定对比
                   number = Inf,            # 提取所有基因
                   sort.by = "P",           # 按 p 值排序
                   adjust.method = "BH")   # FDR 校正方法

# 10. 筛选显著 DEG
deg <- subset(result,
              adj.P.Val < 0.05 & abs(logFC) > 1)
cat("显著差异基因数:", nrow(deg), "\n")
write.csv(result, "limma_results.csv")     # 保存全部结果

批次效应校正（limma 特有优势）¶

# limma 可以直接在设计矩阵中校正批次效应
batch <- factor(c("batch1","batch1","batch2","batch1","batch2","batch2"))

# 在设计矩阵中加入批次作为协变量
design_batch <- model.matrix(~ 0 + group + batch)
colnames(design_batch)[1:2] <- c("HC","T2D")  # 修正列名

# 重新拟合（包含批次校正）
v_batch <- voom(dge, design_batch, plot=FALSE)
fit_batch <- lmFit(v_batch, design_batch)
# ... 后续步骤同上

# 或者用 removeBatchEffect 直接校正表达值（仅用于可视化）
corrected_expr <- removeBatchEffect(v$E, batch = batch, design = design)

参数详解¶

函数/参数	说明	默认值
`voom()`	count → 连续值转换（适合 RNA-seq）	—
`lmFit()`	线性模型拟合	—
`makeContrasts()`	定义组间对比关系	—
`contrasts.fit()`	应用对比矩阵	—
`eBayes()`	经验贝叶斯方差压缩	—
`topTable()`	提取排序结果	number=10
`adjust.method`	FDR 校正方法（BH/BY/none）	BH
`sort.by`	排序依据（P/logFC/AveExpr）	B
`lfc`	最小 log2FC 阈值	0

实战案例：微阵列数据（limma 原始用途）¶

# 微阵列数据直接用 limma（不需要 voom）
# expr_matrix：log2 归一化后的表达矩阵

expr_matrix <- read.table("microarray_norm.txt",
                           header=TRUE, row.names=1, sep="\t")

design_ma <- model.matrix(~ 0 + group)
colnames(design_ma) <- levels(group)

# 直接拟合（不需要 voom，数据已是连续值）
fit_ma <- lmFit(expr_matrix, design_ma)    # 拟合线性模型
fit_ma2 <- contrasts.fit(fit_ma,
    makeContrasts(T2D-HC, levels=design_ma))
fit_ma2 <- eBayes(fit_ma2)                 # 经验贝叶斯

result_ma <- topTable(fit_ma2, n=Inf, adjust="BH")

常见报错与解决¶

报错信息	原因	解决方法
`design matrix not full rank`	设计矩阵存在线性相关	检查并去除冗余列
`voom: 'lib.size' missing`	DGEList 未传入 voom	先创建 DGEList，再传给 voom
`NaN in eBayes`	数据存在极端值	检查并过滤异常样本
`contrasts not estimable`	对比矩阵与设计矩阵不匹配	检查 levels 参数
结果与 DESeq2 差异大	方法不同，两者都合理	取交集作为最可靠结果

速查表¶

# limma-voom 最简流程
library(limma); library(edgeR)
dge <- DGEList(counts, group=group)
dge <- dge[filterByExpr(dge),]
dge <- calcNormFactors(dge)
design <- model.matrix(~0+group)
colnames(design) <- levels(group)
v <- voom(dge, design)
fit <- eBayes(contrasts.fit(lmFit(v,design), makeContrasts(T2D-HC, levels=design)))
res <- topTable(fit, coef="T2D - HC", n=Inf, adjust="BH")
write.csv(res, "limma_results.csv")