缘由
不懂原理的生信人，不是一个好的调包侠。
这个系列不是教程系列 而是原理系列
碎碎念
我的梦想，值得我本人去争取，我今天的生活，绝不是我昨天生活的冷淡抄袭。——《红与黑》
Question
- 如何计算蛋白组中无生物学重复的P值？
- 计算方法的数学原理
Ref
1. MaxQuant 计算蛋白 ratio 的基本方法（重点：中位数）
MaxQuant 手册中写得非常清楚：
Protein ratios are calculated as the median of all SILAC pair ratios that belong to peptides contained in this protein. By using this robust kind of averaging, outlier ratios will leave the protein ratio largely unaffected.
👉 蛋白 ratio = 所有肽段 ratio 的中位数（不是均值）
这是第一层“稳健性”：
中位数不受极端肽段影响。
如果做了归一化（normalization）：
A version of protein ratios that is corrected for unequal total protein amounts is calculated by taking the median of the normalized versions of SILAC pair ratios.
2. 关键问题：没有重复，p 值从哪里来？
The objective in proteomics is usually to determine outlier protein ratios in comparison to the distribution of all protein ratios.
Our scores implement and substantially improve this tradition…
👉 把所有蛋白的 log-ratio 分布看成“零假设下的噪声分布”
👉 如果某个蛋白偏离这个分布足够远，就认为它显著变化
这就是 Significance A/B 能“无生物重复算 p 值”的本质。
🧮 3. 为什么不用传统 Z-score？因为不稳健！
传统 Z-score：
但蛋白质组的 ratio 分布 不是正态分布：
- 有长右尾
- outlier 多
- 低丰度蛋白噪声大
- 高丰度蛋白非常稳定
传统均值 μ 和标准差 σ 会被 outlier拉得面目全非。
所以 MaxQuant 放弃均值/标准差，采用 “分位数” 来构建稳健的 Z-score。
🔧 4. Significance A：用分位数估计“左右标准差”
To make a robust and asymmetrical estimate of the standard deviation…
we calculate the 15.87, 50, and 84.13 percentiles r−1, r0, and r1.
We define r1 − r0 and r0 − r−1 as the right- and left-sided robust standard deviations.
这利用了一个关键性质：
在正态分布中：
- 84.13% = mean + 1 SD
- 50% = mean
- 15.87% = mean − 1 SD
于是 MaxQuant 构造：
右侧稳健 SD：
左侧稳健 SD：
然后：
若 r > r₀（上调）：
若 r < r₀（下调）：
p 值由：
这就是 Significance A。
⚡ 5. Significance A 的问题：忽略了丰度对噪声的影响
- 高丰度蛋白的 ratio 更稳定 → 分布更窄
- 低丰度蛋白噪声大 → 分布更宽
Significance A 用的是“所有蛋白混在一起的分布”，会导致：
- 高丰度蛋白显著性被低估
- 低丰度蛋白显著性被高估
6. Significance B：按丰度分组（binning），分别计算背景分布
英文原文：
Intensity-dependent error is accounted for by calculating the score in bins of similar intensity.
做法：
-
按 intensity / iBAQ 对蛋白排序
-
分成多个 bin（每 bin ≥ 300 个蛋白）
-
在每个 bin 中独立计算
- r₋1、r₀、r₁
- s_+、s_-
- Z
- p
-
最终得到 Significance B（intensity corrected）
本质：
👉 Significance B = 丰度分层后的稳健 Z-score，更准确
7. 可复现 R 代码
Significance A
# ================================
# Significance A (robust Z-score)
# ================================
significance_A <- function(ratio) {
r <- log2(as.numeric(ratio))
qs <- quantile(r, probs = c(0.1587, 0.50, 0.8413), na.rm = TRUE)
rl <- qs[1]
rm <- qs[2]
rh <- qs[3]
sapply(r, function(x){
if (x >= rm) {
z <- (x - rm) / (rh - rm)
} else {
z <- (rm - x) / (rm - rl)
}
pnorm(z, lower.tail = FALSE)
})
}Significance B（丰度分组版）
# ================================
# Significance B (intensity-binned)
# ================================
significance_B <- function(ratio, intensity, min_bin_size = 300) {
df <- data.frame(
ratio = ratio,
log2ratio = log2(as.numeric(ratio)),
intensity = intensity
)
# 按丰度排序
df <- df[order(df$intensity), ]
N <- nrow(df)
nbins <- ceiling(N / min_bin_size)
# 分 bin
df$bin <- cut(seq_len(N), breaks = nbins, labels = FALSE)
# 在每个 bin 中计算 Significance A
df$pvalue_B <- unlist(
tapply(df$ratio, df$bin, function(sub_ratio){
significance_A(sub_ratio)
})
)
return(df)
}总结
① 无生物重复的 p 值来自全体蛋白 ratio 构成的“背景噪声分布”。
② Significance A 用的是全局稳健 Z-score（分位数取代均值/标准差）。
③ Significance B 则按丰度分组，使用“各自的噪声”，更准确。