BioF3 组学数据分析

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

导出日期：2026年6月27日

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 通路激酶变化能直接对应到生物学解释。

报告里要交代的关键参数

常见坑

坑 1：Localization probability 阈值设太松

0.5 的定位概率意味着有一半可能性磷酸化不在这个位点上。下游分析基于错误位点做的 kinase-substrate 映射全是噪声。建议 ≥ 0.75 做探索性分析，≥ 0.9 做最终报告。Class I sites（localization prob ≥ 0.75）是领域共识。

坑 2：没做蛋白丰度归一化就解读磷酸化变化

磷酸位点强度上升可能只是因为蛋白本身丰度增加了，而不是修饰比例改变。不做蛋白丰度归一化（phospho/total protein），会把大量"蛋白量变"误报为"磷酸化调控"。同一组样本必须同时有 total proteome 数据做归一化。

坑 3：KSEA 的激酶-底物数据库覆盖不全

PhosphoSitePlus 对人类主要激酶覆盖还行，但对非模式物种或新发现的激酶几乎没有注释。如果 KSEA 结果只显示 3-4 个激酶有变化，不代表其他激酶没活性变化——可能只是数据库没收录。报告中应注明数据库版本和覆盖率。

坑 4：把 multiplicity > 1 的位点和单磷酸化位点混在一起分析

一个肽段上有两个磷酸化位点（multiplicity = 2）时，MaxQuant 会分别报告。但这两个位点的定量是关联的（来自同一条肽段），不能当作独立观测。通常建议只分析 multiplicity = 1 的位点，除非专门研究多位点共修饰。

下一步

接着深入：

横向延伸：

参考资源