08 与转录组联合分析

蛋白丰度和 mRNA 丰度的相关性其实没有想象中高（典型相关系数 0.3-0.6）。这种不一致本身就是信息：同时有转录组和蛋白组数据时，比较两者能识别"翻译/降解层面被调控"的基因，是单一数据无法看到的故事。

为什么相关性不高

原因	说明
mRNA 半衰期短，蛋白稳定	同一时刻"流量"和"存量"测的是不同东西
翻译效率不同	5'UTR、CDS 长度、密码子偏好
蛋白降解动态	泛素化标记、自噬清除
检测灵敏度差	高丰度蛋白容易测，低丰度的 mRNA 反而灵敏

所以相关性低不是"数据有问题"，是真实的生物学层级。

数据对齐：两份矩阵的统一坐标

转录组的基因 ID 通常是 ENSEMBL，蛋白组多用 UniProt。先映射到 gene symbol 作为公共坐标：

library(biomaRt)
library(dplyr)

mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# RNA: ENSG -> SYMBOL
rna_map <- getBM(c("ensembl_gene_id", "hgnc_symbol"),
                 filters = "ensembl_gene_id",
                 values = rownames(rna_mat),
                 mart = mart)

# Protein: UniProt -> SYMBOL
prot_map <- getBM(c("uniprotswissprot", "hgnc_symbol"),
                  filters = "uniprotswissprot",
                  values = rownames(prot_mat),
                  mart = mart)

# 同时存在的基因
common <- intersect(rna_map$hgnc_symbol, prot_map$hgnc_symbol)
common <- common[common != ""]

相关性图：每个基因一对点

最直观的图：把每个基因的 mRNA logFC 和 protein logFC 画散点。同向变化的落在 y=x 附近，异向的偏离对角线。

library(ggplot2)

joint <- data.frame(
  gene = common,
  rna_logFC = rna_de[common, "logFC"],
  prot_logFC = prot_de[common, "logFC"],
  rna_padj = rna_de[common, "padj"],
  prot_padj = prot_de[common, "padj"]
)

joint$category <- with(joint, case_when(
  rna_padj < 0.05 & prot_padj < 0.05 & sign(rna_logFC) == sign(prot_logFC) ~ "一致变化",
  rna_padj < 0.05 & prot_padj < 0.05 & sign(rna_logFC) != sign(prot_logFC) ~ "反向变化",
  rna_padj < 0.05 & prot_padj >= 0.05 ~ "仅 RNA",
  rna_padj >= 0.05 & prot_padj < 0.05 ~ "仅蛋白",
  TRUE ~ "无显著"
))

ggplot(joint, aes(rna_logFC, prot_logFC, color = category)) +
  geom_point(alpha = 0.6, size = 1) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey60") +
  geom_hline(yintercept = 0, linetype = "dotted", color = "grey60") +
  geom_vline(xintercept = 0, linetype = "dotted", color = "grey60") +
  scale_color_manual(values = c("一致变化" = "#10b981",
                                 "反向变化" = "#f59e0b",
                                 "仅 RNA" = "#2563eb",
                                 "仅蛋白" = "#dc2626",
                                 "无显著" = "grey80")) +
  labs(x = "mRNA log2FC", y = "Protein log2FC") +
  theme_classic()

四种类别各代表不同的生物学：

• 一致变化：转录控制为主
• 反向变化：翻译/降解强烈反向调控（少见但重要）
• 仅 RNA：转录起来了但还没翻译，或蛋白半衰期长来不及响应
• 仅蛋白：稳态层面翻译效率或降解被调控

翻译效率（TE）估计

如果有同样本的核糖体 profiling 数据（Ribo-seq）就能直接算 TE = ribosome footprints / mRNA。没有的话，蛋白丰度 / mRNA 丰度可以作为"稳态翻译效率"的近似：

te <- prot_logFC - rna_logFC  # 注意：log 空间下相减 = 比值

joint$te <- te
ggplot(joint, aes(reorder(gene, te), te)) +
  geom_col(aes(fill = te > 0)) +
  coord_flip() +
  labs(x = "", y = "Protein/RNA ratio change (log2)")

TE 显著上升的基因往往是被 mTOR、eIF 等翻译因子调控的核心节点。

多组学整合方法（跨多种数据类型）

如果再加上代谢组、磷酸化等数据，简单散点不够用。常见整合方法：

方法	适用	R 包
MOFA+	找跨组学共同变化的因子	MOFA2
DIABLO	监督学习，做分类	mixOmics
MultiAssayExperiment	数据结构容器	MultiAssayExperiment

MOFA 是无监督的，给一个隐变量空间，每个 factor 在各组学的 loading 表示"这个因子在 RNA / 蛋白 / 代谢组里分别由哪些 feature 支配"。

library(MOFA2)

mofa <- create_mofa(list(RNA = rna_mat, Protein = prot_mat))
mofa <- prepare_mofa(mofa,
  data_options = list(scale_views = TRUE),
  model_options = list(num_factors = 10),
  training_options = list(seed = 42)
)
mofa <- run_mofa(mofa, save_data = FALSE)

plot_variance_explained(mofa)
plot_factor(mofa, factors = 1:3)

报告里要交代

• 数据来自同样本还是不同样本（影响相关性解读）
• ID 映射的版本（biomaRt / ENSEMBL release 号）
• 一致变化/仅蛋白/仅 RNA 各有多少
• 翻译效率改变的 top 基因

参考资源

• Liu 等的综述："On the dependency of cellular protein levels on mRNA abundance"（Cell 2016）
• MOFA2 文档：biofam.github.io/MOFA2
• mixOmics 教程：DIABLO 与多块整合