BioF3 组学数据分析

08 与转录组联合分析

导出日期：2026年6月27日

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()

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

翻译效率（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)

报告里要交代

常见坑

坑 1：看到 RNA-protein 不一致就说"翻译调控"

相关性低的原因很多：时间延迟（mRNA 先变，蛋白 12-48 小时后才跟上）、蛋白半衰期长、检测灵敏度差异。不能看到散点图偏离对角线就下"翻译后调控"的结论。需要排除时间因素（同一时间点的样本）和技术因素（蛋白检测灵敏度）后才能讨论生物学机制。

坑 2：ID 映射版本不一致导致假"不一致"

RNA-seq 用 Ensembl release 105 的 gene symbol，蛋白组用 UniProt 2024 版映射。两个版本对同一个蛋白的 symbol 可能不同（比如历史上改过名的基因）。这会让本来一致的变化因为 ID 对不上而被误判为"仅 RNA"或"仅蛋白"。

坑 3：只看单基因 logFC 散点忽略了统计不确定性

RNA logFC = 2（padj = 0.001）和 protein logFC = 0.5（padj = 0.3），画在散点图上看起来"不一致"，但蛋白的 logFC 置信区间可能包含 2。不显著不代表真的没变化，可能只是检测力不够。建议用颜色区分"双层显著""单层显著""都不显著"。

坑 4：翻译效率（TE）的近似估计被当作精确值

protein logFC - RNA logFC 作为 TE 变化的近似，假设蛋白降解速率不变且两层定量误差独立。但蛋白降解是动态的，这个近似在稳态改变剧烈时（如应激响应）会失效。只有 Ribo-seq 数据才能真正量化翻译效率。

下一步

接着深入：

横向延伸：

参考资源