很多人做转录组（RNA-seq）分析，跑完差异表达基因、拉个火山图、做做GO/KEGG富集，就觉得完工了。但这时候你手里拿到的，只是一堆“名字”，至于这些基因到底在细胞里干了什么、怎么调控的、上游信号是谁、下游效应是什么，还是懵的。
这时候，蛋白质组学（Proteomics）和转录组学的联合分析，就成了破局的关键。
今天不聊那些高大上的论文术语，咱们就把这件事当成一个“从线索到真相”的侦探故事来讲。我会带你一步步拆解：为什么RNA数据不够用？怎么把RNA和蛋白数据配对？常见坑有哪些？最后给你一个完整的解读思路，让你下次跑完多组学数据，能自信地说：“这波分析，我懂。”
一、先别急，为什么RNA够用了，我还要蛋白？
咱们先来打个比方。
想象你在一家餐厅后厨观察。RNA就像厨房里的订单小票——上面写着“要做10份宫保鸡丁”“5份红烧肉”。订单多，说明老板觉得这道菜要火，准备多备料。
但订单多，不等于菜就一定做出来了，更不代表端上桌的味道好不好，更不代表顾客吃不吃得完。
- 可能厨师没空做（翻译抑制）
- 可能做完又扔了（蛋白降解）
- 可能做好了，但顾客不吃（功能未执行）
这就是中心法则的“断层”：mRNA水平 ≠ 蛋白水平。
举个真实例子
假设你研究癌症药物耐药性。RNA-seq显示，某个凋亡基因BAX的mRNA上调了2倍，你以为是“细胞要凋亡了，药有效”。
但蛋白水平检测发现，BAX蛋白反而下降了3倍。
为什么？因为细胞启动了反馈机制，通过泛素-蛋白酶体系统快速降解了这个蛋白，抵抗死亡信号。
如果你只看RNA，结论完全反了。
所以，RNA告诉你“计划做什么”，蛋白告诉你“实际做了什么”。两者结合，才能还原真实的生物学状态。
二、蛋白质组学数据从哪来？三种主流技术
在开始分析之前，你得知道你的蛋白数据是怎么来的。目前主流有三种：
1. 质谱蛋白质组学（LC-MS/MS）—— 无偏发现
这是目前最主流的方法。原理是把蛋白切 peptides，通过质谱仪测出肽段的质荷比，再反推蛋白身份和丰度。
- 优点：无偏，能发现新蛋白、修饰位点、剪接变体
- 缺点：动态范围有限，低丰度蛋白难检测；定量精度依赖同位素标记（如TMT、SILAC）或无标记（Label-free）
- 适合：探索性研究，发现差异蛋白
2. 靶向蛋白组学（SRM/MRM、PRM）—— 精准验证
在发现阶段找到候选蛋白后，用这种方法精准定量特定几个蛋白。
- 优点：灵敏度极高，定量准确，可重复性好
- 缺点：只能检测预设的蛋白，不能发现新东西
- 适合：临床验证、生物标志物确认
3. 抗体芯片/MSD/Ella —— 高通量但有限
- 优点：操作简单，通量高，适合大量样本
- 缺点：依赖抗体质量，只能测已知蛋白，覆盖有限
- 适合：临床样本筛选、机制验证
建议：如果你是做科研探索，首选LC-MS/MS；如果是要发高分文章，最好配上靶向验证或功能实验。
三、RNA+蛋白联合分析的核心流程
别被流程图吓到，其实逻辑很简单：配对、校正、整合、解读。
第一步：数据预处理与质量控制
RNA和蛋白数据来自不同批次、不同平台，必须先做质控。
RNA-seq质控要点
- 检查测序深度：通常建议>30M clean reads
- 检查重复样本相关性：Pearson相关系数应>0.9
- 检查PCA图：看是否有批次效应，样本是否按分组聚类
蛋白质组学质控要点
- 检查缺失值比例：如果某个蛋白在>50%样本中缺失，考虑剔除
- 检查技术重复：CV值应<20%
- 检查峰面积分布：看是否有异常高点（可能是污染）
实操建议：用R语言画几个图，心里有底再往下走。
# RNA-seq样本相关性矩阵热图
library(pheatmap)
cor_matrix <- cor(log2(counts + 1), method = "pearson")
pheatmap(cor_matrix, clustering_distance_rows = "correlation",
clustering_method = "ward.D2", show_rownames = TRUE,
color = colorRampPalette(c("blue", "white", "red"))(100))
# 蛋白数据缺失值分布
missing_rate <- colMeans(is.na(protein_data))
barplot(missing_rate, main = "Missing Value Rate per Protein",
ylim = c(0, 1), ylab = "Missing Rate", las = 2)
abline(h = 0.5, col = "red", lty = 2)
第二步：数据标准化与批次校正
RNA和蛋白的数据分布不一样，必须标准化。
RNA标准化
- 常用方法：TPM、FPKM、或DESeq2的median ratio法
- 注意：不要直接用原始count做下游分析
蛋白标准化
- 常用方法：-log2转换 + median中心化
- 批次效应：用ComBat或limma的removeBatchEffect
RNA-蛋白联合标准化难点
RNA和蛋白的量纲完全不同，不能直接比较数值。通常的做法是：
- 分别对RNA和蛋白数据做Z-score标准化（使均值为0，标准差为1）
- 或者分别做rank归一化，然后比较rank相关性
# Z-score标准化示例
rna_z <- t(scale(t(rna_data)))
protein_z <- t(scale(t(protein_data)))
# 计算RNA与蛋白的相关性
cor_matrix <- cor(rna_z, protein_z, method = "pearson")
heatmap(cor_matrix, col = cm.colors(256),
main = "RNA-Protein Correlation Heatmap")
第三步：差异分析与整合
分别对RNA和蛋白做差异分析，找出各自显著变化的分子。
RNA差异分析（DESeq2示例）
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = count_table,
colData = sample_info,
design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "treated", "control"))
res_sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
蛋白差异分析（limma示例）
library(limma)
design <- model.matrix(~ 0 + condition)
colnames(design) <- c("Control", "Treated")
fit <- lmFit(protein_z, design)
contrast.matrix <- makeContrasts(Treated - Control, levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
protein_sig <- topTable(fit2, number = Inf, adjust.method = "BH")
protein_sig <- subset(protein_sig, adj.P.Val < 0.05 & abs(logFC) > 1)
整合策略
- 一致变化：RNA↑ + Protein↑（最可信）
- 转录后调控：RNA↑ + Protein↓（或反之）（最值得深挖）
- 单向变化：只有RNA或只有蛋白变化（需验证）
# 合并结果
combined <- merge(rna_sig, protein_sig,
by = "gene", suffixes = c("_RNA", "_Protein"),
all = TRUE)
# 分类
combined$type <- "Consistent_Up"
combined$type[is.na(combined$RNA_log2FC) & !is.na(combined$Protein_log2FC)] <- "Protein_only"
combined$type[!is.na(combined$RNA_log2FC) & is.na(combined$Protein_log2FC)] <- "RNA_only"
combined$type[combined$RNA_log2FC > 0 & combined$Protein_log2FC < 0] <- "Post_transcriptional"
combined$type[combined$RNA_log2FC < 0 & combined$Protein_log2FC > 0] <- "Post_transcriptional"
# 可视化
library(ggplot2)
ggplot(combined, aes(x = RNA_log2FC, y = Protein_log2FC, color = type)) +
geom_point(alpha = 0.6) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
geom_hline(yintercept = c(-1, 1), linetype = "dashed") +
theme_minimal() +
labs(title = "RNA vs Protein Expression Integration",
x = "RNA log2FoldChange", y = "Protein log2FoldChange")
第四步：功能注释与通路分析
找到差异分子后，下一步是“这些分子在干嘛”。
常用工具
- GO分析：Gene Ontology，描述基因产物参与的生物过程、分子功能、细胞组分
- KEGG通路：代谢通路、信号通路
- Reactome：更详细的反应级联
- STRING：蛋白互作网络
- GSEA：基因集富集分析，看整个通路是否协同变化
RNA+蛋白联合富集的亮点
传统GO分析只分析RNA差异基因，联合分析可以：
- 提高可信度：RNA+蛋白一致变化的基因，功能更可靠
- 发现新机制：转录后调控相关的通路（如泛素化、自噬）会被凸显
- 网络视角：构建RNA-蛋白共调控网络，找到枢纽分子
# 使用clusterProfiler做GO分析
library(clusterProfiler)
library(org.Hs.eg.db)
# RNA差异基因GO
ego_rna <- enrichGO(gene = sig_rna_genes,
OrgDb = org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05)
# 蛋白差异基因GO
ego_prot <- enrichGO(gene = sig_protein_genes,
OrgDb = org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05)
# 比较两个结果
barplot(ego_rna, showCategory = 20) + ggtitle("RNA GO Terms")
barplot(ego_prot, showCategory = 20) + ggtitle("Protein GO Terms")
四、深度解读：如何从“数据”到“故事”
做完分析只是第一步，真正考验功力的是解读。
1. 找枢纽分子
在蛋白互作网络（PPI）中，那些连接度最高的节点，往往是关键调控因子。
# 使用STRING数据库构建PPI网络
library(STRINGdb)
string <- STRINGdb$new(version = "12.0", species = 9606, score_threshold = 400)
genes <- c("TP53", "AKT1", "MYC", "EGFR", "VEGFA", "STAT3") # 示例基因
net <- string$get_network(genes)
plot(net)
2. 追踪信号通路
选一条你觉得重要的通路，比如MAPK或PI3K/AKT通路，把通路上所有分子的RNA和蛋白数据拉出来，画成通路图。
你可以用Pathway Commons、Reactome或者Ingenuity Pathway Analysis (IPA)来做。
关键问题：
- 哪些节点是RNA变化但蛋白没变？可能是转录后调控
- 哪些节点是蛋白变化但RNA没变？可能是翻译调控或蛋白降解
- 上游信号是谁？下游效应是什么？
3. 跨组学一致性评估
计算RNA与蛋白的相关系数，看整体一致性有多高。
- 高一致性（r > 0.7）：转录主导，蛋白水平主要受mRNA调控
- 中等一致性（r = 0.4-0.7）：转录+转录后双重调控
- 低一致性（r < 0.4）：翻译后调控主导（如蛋白修饰、降解）
这个分析本身就能回答很多机制问题。
五、常见坑与解决方案
坑1：样本量太小
问题：蛋白组学成本远高于RNA-seq，很多人只测3-5个生物学重复。
后果：统计效力低，假阳性/假阴性率高。
解决：
- 尽量保证≥6个生物学重复
- 使用paired设计（如果可能）
- 采用更严格的阈值（padj < 0.01，log2FC > 1.5）
坑2：缺失值处理不当
问题：蛋白数据中缺失值很多，直接删除会导致信息丢失。
解决：
- 使用kNN或minProb等方法填补低丰度缺失值
- 对于系统性缺失（只在某些样本中检测不到），考虑删除
- 报告缺失值比例，让审稿人知道数据质量
坑3：RNA与蛋白匹配错误
问题：RNA和蛋白数据来自不同批次、不同平台，甚至不同个体。
解决：
- 尽量使用同一批样本分别测RNA和蛋白
- 如果无法匹配，用paired statistical methods（如MOFA+）
- 在方法部分详细说明样本来源和匹配方式
坑4：忽略翻译后修饰
问题：蛋白功能往往由修饰决定（磷酸化、乙酰化等），但普通蛋白组学只测总蛋白。
解决：
- 如果预算允许，加做磷酸化蛋白组学（Phosphoproteomics）
- 利用现有数据中的修饰位点信息（如UniProt注释）
- 结合RNA-seq和蛋白组学推断调控机制
坑5：过度解读低置信度结果
问题：RNA和蛋白都显著变化，就认为是“关键靶点”。
解决：
- 验证：qPCR验证RNA，Western blot验证蛋白
- 功能实验：敲低/过表达，看表型是否改变
- 文献支持：查查这个分子在同类研究中是否有报道
六、实战案例：从数据到结论
让我讲一个完整的分析故事，帮你把零散的知识点串起来。
研究背景
研究某种新型抗癌药物对乳腺癌细胞系的药效。
实验设计
- 分组：对照组 vs 药物处理组
- 重复：6个生物学重复
- 技术：RNA-seq + Label-free蛋白组学
- 时间点：处理24小时
分析流程
1. 质控
- RNA-seq：平均测序深度45M reads，RiboAlert评分良好，PCA显示组内聚集，组间分离
- 蛋白组：平均覆盖蛋白5000个，技术重复CV=12%，缺失值<30%
2. 差异分析
- RNA：找到850个差异表达基因（DEGs）
- 蛋白：找到420个差异蛋白（DEPs）
3. 整合分析
- 一致上调：350个基因/蛋白
- 一致下调：180个基因/蛋白
- 转录后调控（RNA↑蛋白↓）：70个
- 转录后调控（RNA↓蛋白↑）：45个
4. 功能富集
- 一致变化基因显著富集在“细胞周期”、“DNA复制”、“有丝分裂”
- 转录后调控基因显著富集在“泛素介导的蛋白水解”、“蛋白修饰”
- KEGG分析显示“p53 signaling pathway”、“cell cycle”显著激活
5. 网络分析
- PPI网络识别出CDK1、CCNB1、PLK1为枢纽节点
- 这些节点在RNA和蛋白水平均显著上调，且互作强度高
6. 验证
- qPCR验证：随机选10个基因，验证结果与测序数据一致
- Western blot：验证关键蛋白CDK1、CCNB1表达上调
- 细胞实验：药物处理导致细胞周期G2/M期阻滞，与通路分析一致
结论
该药物通过激活p53通路，上调细胞周期调控蛋白，导致乳腺癌细胞G2/M期阻滞，发挥抗癌作用。
七、工具推荐与资源汇总
分析工具
| 用途 | 推荐工具 |
|---|---|
| RNA-seq分析 | DESeq2, edgeR, limma-voom |
| 蛋白组分析 | MaxQuant, DIA-NN, Proteome Discoverer |
| 差异蛋白分析 | limma, msGEO, DEP |
| 整合分析 | mixOmics, MOFA+, Paintomics |
| 功能富集 | clusterProfiler, GSEA, Enrichr |
| 网络分析 | Cytoscape, STRING, gephi |
| 可视化 | ggplot2, pheatmap, ComplexHeatmap |
数据库
- STRING：蛋白互作网络
- KEGG：代谢与信号通路
- Reactome：详细生物反应
- UniProt：蛋白注释与修饰信息
- CPTAC：癌症蛋白质组学数据库
- PRIDE：蛋白质组学原始数据仓库
八、给你的几点建议
- 实验设计先行：在测序之前就想好怎么整合分析