06 翻译后修饰：磷酸化位点定量

蛋白量变只能告诉你"多了"或"少了"，但很多生物学过程靠的是修饰：磷酸化、乙酰化、泛素化等。质谱在富集后能定量到位点级别，本章以最常见的磷酸化为例。

富集 + 检测流程

普通 shotgun 蛋白质组识别不到磷酸肽。需要先做磷酸肽富集（IMAC、TiO₂、Fe-NTA），再上质谱：

步骤	工具	输出
富集 + 上机	TiO₂、IMAC	磷酸化富集的肽段
数据库搜索	MaxQuant (PTM) / FragPipe	Phospho (STY)Sites.txt
位点定位	PTM Score / AScore	每个位点定位概率
定量	LFQ / TMT	位点强度矩阵
差异分析	limma / DEqMS	差异磷酸位点
通路解释	KSEA / WebGestalt	上游激酶活性变化

起点：Phospho (STY)Sites.txt

MaxQuant 输出的位点表格每一行是一个位点（不是蛋白），列里有 Position、Sequence window、Localization prob 等。常用过滤：

library(dplyr)
library(readr)

ph <- read_tsv("Phospho (STY)Sites.txt") %>%
  filter(`Localization prob` >= 0.75,        # 位点定位概率
         is.na(Reverse),                      # 去除反向库
         is.na(`Potential contaminant`),      # 去除污染
         `Localization prob` >= 0.75)

# 提取定量列（multiplicity = 1 的位点）
intensity_cols <- grep("^Intensity ", names(ph), value = TRUE)
mat <- ph %>% select(all_of(intensity_cols)) %>% as.matrix()
rownames(mat) <- paste0(ph$`Gene names`, "_", ph$`Amino acid`, ph$Position)

定位概率 0.75 是常用阈值，更严格可以 0.9。

蛋白丰度归一化

磷酸位点的强度变化可能来自蛋白本身丰度变化（蛋白多了，磷酸位点自然多），也可能是修饰比例变化。要分清楚，需要在同样本上同时测蛋白丰度（proteinGroups.txt），然后归一化：

# ph_mat: 位点强度矩阵
# prot_mat: 同一组样本的蛋白丰度矩阵（log2 强度）
# protein_id: ph 表中每个位点对应的蛋白 ID

ph_norm <- ph_mat
for (i in seq_len(nrow(ph_mat))) {
  pid <- protein_id[i]
  if (pid %in% rownames(prot_mat)) {
    ph_norm[i, ] <- ph_mat[i, ] - prot_mat[pid, ]
  }
}

归一化后的 ph_norm 反映的是"磷酸化比例"的变化，不再受蛋白丰度干扰。

差异分析：limma + 位点权重

DEqMS 在 limma 基础上加入位点强度权重，对低强度位点更稳健：

library(limma)
library(DEqMS)

design <- model.matrix(~0 + group, data = sample_info)
colnames(design) <- levels(sample_info$group)
contrast <- makeContrasts(Treatment - Control, levels = design)

fit <- lmFit(ph_norm, design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)

# DEqMS 加位点 PSM 数权重
fit2$count <- ph$`Number of PSM` # 每个位点的支持 PSM 数
fit_deqms <- spectraCounteBayes(fit2)

result <- outputResult(fit_deqms, coef_col = 1)

上游激酶富集（KSEA / KSEAapp）

差异位点列表本身意义有限，更有用的是问"哪些激酶被激活了"。KSEA 把每个磷酸位点映射到已知激酶 → 底物关系（来自 PhosphoSitePlus），再做 GSEA-like 富集：

library(KSEAapp)

ksea_input <- result %>%
  transmute(Protein = gene,
            Gene = gene,
            Peptide = sequence,
            Residue.Both = paste0(amino_acid, position),
            p = adj.P.Val,
            FC = logFC)

ksea <- KSEA.Complete(ksea_input,
                       NetworKIN = TRUE,
                       NetworKIN.cutoff = 5,
                       m.cutoff = 5,
                       p.cutoff = 0.05)

输出每个激酶的活性 z-score 和 p 值。常见的 PI3K-AKT、MAPK、CDK 通路激酶变化能直接对应到生物学解释。

报告里要交代的关键参数

• 富集方法（IMAC / TiO₂）和富集次数
• Localization probability 阈值（推荐 0.75）
• 是否做了蛋白丰度归一化
• 激酶 → 底物数据库版本（PhosphoSitePlus 更新很频繁）

参考资源

• MaxQuant PTM workflow 官方文档
• DEqMS 包文档：质谱定量数据的 limma 扩展
• PhosphoSitePlus 数据库：phosphosite.org
• KSEAapp R 包：激酶活性富集分析