蛋白质组差异分析limma¶

传统t检验的问题（小样本时）：
  - 每组3个重复 → 方差估计不准
  - 一个蛋白质方差很小 → 假阳性
  - 一个蛋白质方差很大 → 假阴性

limma的解决办法（经验贝叶斯）：
  "如果你一个人估计不准，就看看大家"

  ├── 每个蛋白质有自己的方差估计
  ├── 所有蛋白质的方差形成一个"大家庭"
  ├── 用"大家庭"的平均水平来调整每个蛋白质的方差
  └── 结果：小样本也能得到稳定的差异分析结果

  打比方：
  你班上只有3个人考试，很难判断平均水平
  但如果参考全年级1000人的成绩，就能更准确地判断你们班的水平

#!/usr/bin/env Rscript
# 蛋白质组学limma差异分析完整流程

# ========== 安装和加载包 ==========
if (!requireNamespace("BiocManager", quietly = TRUE))  # 检查BiocManager
    install.packages("BiocManager")
BiocManager::install("limma")  # 安装limma

library(limma)        # 差异分析
library(ggplot2)      # 画图
library(dplyr)        # 数据处理

# ========== 第一步：读取MaxQuant结果 ==========
pg <- read.delim(
    "combined/txt/proteinGroups.txt",  # MaxQuant蛋白质组输出
    stringsAsFactors = FALSE  # 字符串不自动转因子
)
cat("原始蛋白质组数:", nrow(pg), "\n")  # 打印原始数量

# ========== 第二步：数据过滤 ==========
# 去除反库、污染物、仅位点鉴定
pg <- pg[pg$Reverse != "+", ]         # 去反库匹配
pg <- pg[pg$Potential.contaminant != "+", ]  # 去污染物
pg <- pg[pg$Only.identified.by.site != "+", ]  # 去仅位点鉴定
cat("过滤后蛋白质数:", nrow(pg), "\n")

# ========== 第三步：提取LFQ定量矩阵 ==========
# 提取LFQ intensity列
lfq_cols <- grep("^LFQ.intensity\\.", colnames(pg), value = TRUE)  # 找LFQ列
lfq_matrix <- as.matrix(pg[, lfq_cols])  # 转为矩阵
rownames(lfq_matrix) <- pg$Gene.names  # 用基因名做行名

# 0替换为NA
lfq_matrix[lfq_matrix == 0] <- NA  # MaxQuant中0表示未检测到

# Log2转换
lfq_log2 <- log2(lfq_matrix)  # 对LFQ值取log2

cat("样品数:", ncol(lfq_log2), "\n")
cat("蛋白质数:", nrow(lfq_log2), "\n")

# ========== 第四步：缺失值过滤 ==========
# 至少在一个组中有70%以上有效值
# 假设前3列是Control，后3列是Treatment
n_ctrl <- 3   # 对照组样品数
n_treat <- 3  # 处理组样品数

# 计算每个蛋白质在每组中的有效值比例
ctrl_valid <- rowSums(!is.na(lfq_log2[, 1:n_ctrl])) / n_ctrl  # 对照组有效率
treat_valid <- rowSums(!is.na(lfq_log2[, (n_ctrl+1):(n_ctrl+n_treat)])) / n_treat  # 处理组有效率

keep <- ctrl_valid >= 0.7 | treat_valid >= 0.7  # 至少一组70%有效
lfq_filtered <- lfq_log2[keep, ]  # 过滤
cat("缺失值过滤后:", nrow(lfq_filtered), "个蛋白质\n")

# ========== 第五步：缺失值填充 ==========
# 用最小值向下偏移的方法填充（模拟低丰度蛋白）
impute_min_shift <- function(mat, shift = 1.8, scale = 0.3) {
    # shift: 从均值向下偏移的标准差数
    # scale: 填充分布的宽度
    global_mean <- mean(mat, na.rm = TRUE)  # 全局均值
    global_sd <- sd(mat, na.rm = TRUE)  # 全局标准差

    fill_mean <- global_mean - shift * global_sd  # 填充均值
    fill_sd <- scale * global_sd  # 填充标准差

    for (i in seq_len(ncol(mat))) {  # 遍历每列
        na_idx <- is.na(mat[, i])  # 找缺失值
        n_na <- sum(na_idx)  # 缺失值数量
        if (n_na > 0) {
            mat[na_idx, i] <- rnorm(n_na, fill_mean, fill_sd)  # 用正态分布随机填充
        }
    }
    return(mat)
}

set.seed(42)  # 设置随机种子（保证可重复性）
lfq_imputed <- impute_min_shift(lfq_filtered)  # 填充缺失值
cat("缺失值填充完成\n")

# ========== 第六步：limma差异分析 ==========
# 创建分组因子
group <- factor(c(rep("Control", n_ctrl), rep("Treatment", n_treat)))  # 定义分组
cat("分组:", levels(group), "\n")

# 创建设计矩阵（告诉limma样品分组信息）
design <- model.matrix(~ 0 + group)  # 无截距模型
colnames(design) <- levels(group)  # 列名=组名
cat("设计矩阵:\n")
print(design)

# 创建对比矩阵（告诉limma比较哪两组）
contrast <- makeContrasts(
    Treatment_vs_Control = Treatment - Control,  # 处理组 vs 对照组
    levels = design
)
cat("对比矩阵:\n")
print(contrast)

# 拟合线性模型
fit <- lmFit(lfq_imputed, design)  # 第一步：拟合线性模型
fit2 <- contrasts.fit(fit, contrast)  # 第二步：应用对比
fit3 <- eBayes(fit2)  # 第三步：经验贝叶斯调节（关键！）

# 提取差异分析结果
results <- topTable(
    fit3,
    coef = "Treatment_vs_Control",  # 比较
    number = Inf,  # 返回所有蛋白质
    sort.by = "P"  # 按p值排序
)

cat("\n===== 差异分析结果摘要 =====\n")
cat("上调蛋白质(log2FC>1, padj<0.05):", 
    sum(results$logFC > 1 & results$adj.P.Val < 0.05), "\n")
cat("下调蛋白质(log2FC<-1, padj<0.05):", 
    sum(results$logFC < -1 & results$adj.P.Val < 0.05), "\n")

# 保存结果
write.csv(results, "limma_DEP_results.csv", row.names = TRUE)  # 保存为CSV
cat("结果已保存到 limma_DEP_results.csv\n")

# ========== 第七步：火山图 ==========
results$significance <- "Not Sig"  # 默认不显著
results$significance[results$logFC > 1 & results$adj.P.Val < 0.05] <- "Up"  # 上调
results$significance[results$logFC < -1 & results$adj.P.Val < 0.05] <- "Down"  # 下调

p <- ggplot(results, aes(x = logFC, y = -log10(adj.P.Val), color = significance)) +
    geom_point(alpha = 0.6, size = 1.5) +  # 散点图
    scale_color_manual(values = c("Up" = "red", "Down" = "blue", "Not Sig" = "gray70")) +  # 颜色
    geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray50") +  # p值阈值线
    geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "gray50") +  # FC阈值线
    labs(x = "Log2 Fold Change", y = "-Log10 Adjusted P-value",
         title = "Volcano Plot - Limma DEP Analysis") +  # 标签
    theme_classic() +  # 经典主题
    theme(legend.position = "top")  # 图例在上方

ggsave("volcano_limma.png", p, width = 8, height = 7, dpi = 300)  # 保存图片
cat("火山图已保存到 volcano_limma.png\n")

# ========== 多组比较（>2组） ==========

# 假设有3个组
group <- factor(c("WT", "WT", "WT", "KO", "KO", "KO", "Rescue", "Rescue", "Rescue"))

# 设计矩阵
design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)

# 多个对比
contrast <- makeContrasts(
    KO_vs_WT = KO - WT,           # KO vs 野生型
    Rescue_vs_KO = Rescue - KO,    # 拯救 vs KO
    Rescue_vs_WT = Rescue - WT,    # 拯救 vs 野生型
    levels = design
)

# 依次分析
fit <- lmFit(lfq_imputed, design)
fit2 <- contrasts.fit(fit, contrast)
fit3 <- eBayes(fit2)

# 提取各比较的结果
res_ko_wt <- topTable(fit3, coef = "KO_vs_WT", number = Inf)      # KO vs WT
res_res_ko <- topTable(fit3, coef = "Rescue_vs_KO", number = Inf)  # Rescue vs KO
res_res_wt <- topTable(fit3, coef = "Rescue_vs_WT", number = Inf)  # Rescue vs WT

# F检验：哪些蛋白质在任意比较中显著
res_ftest <- topTable(fit3, number = Inf)  # 综合F检验
cat("任意比较中显著的蛋白质:", sum(res_ftest$adj.P.Val < 0.05), "\n")

# ========== 差异蛋白质热图 ==========
library(pheatmap)  # 热图包

# 提取显著差异蛋白质
sig_proteins <- rownames(results)[results$adj.P.Val < 0.05 & abs(results$logFC) > 1]
sig_matrix <- lfq_imputed[sig_proteins, ]  # 提取显著蛋白的表达矩阵

# Z-score标准化（按行）
sig_zscore <- t(scale(t(sig_matrix)))  # 对每个蛋白质标准化

# 注释信息
anno_col <- data.frame(
    Group = group,  # 样品分组
    row.names = colnames(sig_zscore)  # 样品名
)

# 画热图
pheatmap(
    sig_zscore,
    annotation_col = anno_col,  # 列注释
    color = colorRampPalette(c("blue", "white", "red"))(100),  # 蓝白红色板
    clustering_method = "ward.D2",  # 聚类方法
    show_rownames = nrow(sig_zscore) < 50,  # 蛋白质太多就不显示行名
    main = "Differentially Expressed Proteins",  # 标题
    filename = "heatmap_DEP.png",  # 保存文件
    width = 8, height = 10  # 图片大小
)
cat("热图已保存到 heatmap_DEP.png\n")

========================================
蛋白质组差异分析limma 速查表
========================================

【limma核心三步】
1. lmFit()           → 拟合线性模型
2. contrasts.fit()   → 应用对比（比较哪两组）
3. eBayes()          → 经验贝叶斯调节（核心！）

【结果提取】
topTable()           → 获取差异分析结果
  coef               → 指定哪个对比
  number = Inf       → 返回所有蛋白质
  sort.by = "P"      → 按p值排序

【结果列含义】
logFC                → log2倍数变化
AveExpr              → 平均表达量
t                    → t统计量
P.Value              → 原始p值
adj.P.Val            → 校正后p值（BH方法）
B                    → B统计量（log-odds）

【常用阈值】
显著差异             → adj.P.Val < 0.05 & |logFC| > 1
强差异               → adj.P.Val < 0.01 & |logFC| > 2

【数据预处理】
过滤                 → 去反库、污染物、仅位点鉴定
Log2转换             → log2(LFQ_intensity)
缺失值过滤           → 至少一组70%有效值
缺失值填充           → 最小值下移法 or KNN
标准化               → 中位数中心化 or quantile

【面试考点】
Q: limma的经验贝叶斯是什么？
A: 用所有蛋白质的方差信息来调节单个蛋白质的方差估计

Q: 为什么不直接用t检验？
A: 小样本时方差不稳定，limma的贝叶斯方法更稳健

Q: 蛋白质组学为什么要做log2转换？
A: 原始LFQ值是右偏分布，log2后接近正态分布，适合参数统计
========================================