正文

蛋白质转录组分析实战从0到1掌握RNA测序与蛋白表达关联数据处理流程及常见坑点避坑指南

/0 浏览量
0905

嘿，朋友！我是你的老朋友 Agnes。今天咱们不聊虚的，直接切入生物信息学里最让人“又爱又恨”的交叉领域——转录组与蛋白质组的联合分析。

我知道你在哪类坑里摔过跟头：明明RNA-seq数据跑得很漂亮，PCA分离得清清楚楚，结果一拉到和质谱（Mass Spectrometry, MS）数据对比的时候，相关系数低得让人怀疑人生；或者是在做差异表达分析时，发现“转录本涨了，蛋白却没涨”，一脸懵圈。

别急，这不是你技术不行，而是中心法则在真核生物里从来都不是简单的 \(DNA \rightarrow RNA \rightarrow Protein\) 直线推导。这里有剪接、降解、翻译效率、蛋白稳定性……无数的“中间商”在赚差价。

这篇指南，我将带你从最基础的数据预处理开始，一步步走到最终的整合分析，并且我会把你可能踩过的每一个坑都挖出来，涂上显眼的警示牌。咱们用大白话讲清楚那些复杂的生物统计学逻辑。

第一部分：为什么RNA和蛋白总是“貌合神离”？

在动手写代码之前，你必须先在脑海里建立这个认知：RNA水平并不等于蛋白水平。

很多初学者（甚至部分老手）会直接拿FPKM或TPM值和蛋白丰度去算Pearson相关系数，发现 \(R^2\) 只有0.4-0.6，然后就觉得数据有问题。其实，这在大多数生物系统中是正常的。

根据经典的蛋白质组学研究（如Cox et al., Science 2014; Bateman et al., Mol Syst Biol 2016），全基因组范围内，mRNA与蛋白的表达相关性中位数大约在 0.4 到 0.7 之间。这意味着什么？意味着只有40%-70%的蛋白变异可以由mRNA变异解释，剩下的30%-60%归因于：

  1. 翻译调控：microRNA、RNA结合蛋白、mRNA二级结构影响核糖体加载效率。
  2. 蛋白降解：泛素-蛋白酶体系统、自噬等途径决定蛋白半衰期。
  3. 翻译后修饰：磷酸化、糖基化等不改变蛋白总量但改变其功能状态。
  4. 技术噪音：RNA-seq和质谱（LC-MS/MS）本身的检测动态范围和误差来源不同。

避坑指南 #1：不要追求接近0.9的RNA-蛋白相关性，除非你在研究一种极其特殊的、翻译偶联极紧密的系统（如酵母在快速生长状态下）。如果相关性低，先检查是生物学现象还是技术批次效应，而不是盲目剔除数据。

第二部分：数据源——你手里都有哪些“武器”？

在进行关联分析前，我们要明确两种数据的技术特性，因为它们的“语言”完全不同。

2.1 RNA-seq 数据的特点

  • 单位：通常是 Counts（原始读数），经过标准化后变成 FPKM/RPKM（旧式，不推荐跨样本比较）或 TPM（Transcripts Per Million，推荐）。
  • 优势：灵敏度极高，能检测到低丰度转录本，能发现新异构体、融合基因、SNP等。
  • 劣势：是间接测量，且受文库制备、测序深度影响大。
  • 关键点：必须做基因水平的整合，不能拿ISOFORM（异构体）直接去配蛋白，除非你的质谱数据能特异区分肽段对应的异构体（这很难）。

2.2 蛋白质组数据（LC-MS/MS）的特点

  • 单位：通常是 LFQ Intensity（Label-Free Quantification，无标记定量强度）或 iBAQ（Intensity Based Absolute Quantification）。iBAQ理论上是绝对定量，更适合跨样本比较绝对丰度。
  • 优势：直接反映功能执行者，更接近表型。
  • 劣势：动态范围比RNA-seq窄（大约3-4个数量级 vs RNA的6-8个），低丰度蛋白难检测，存在“漏值”（Missing Values）问题。
  • 关键点：质谱数据中，一个蛋白往往对应多条肽段，需要通过算法（如MaxQuant）汇总成一个蛋白级别的表达值。

避坑指南 #2：单位对齐。千万不要拿RNA的TPM直接和质谱的LFQ Intensity做图！两者分布形态完全不同（TPM近似对数正态，LFQ往往偏态更严重）。在整合前，两者都需要做Log2转换。

第三部分：实战流程——从原始数据到关联图谱

假设你已经有了两批数据：

  • sample_RNA.bam 或 sample_RNA.fastq
  • sample_protein.mgf 或 MaxQuant输出的 proteinGroups.txt

我们将流程分为五个核心步骤。

步骤一：RNA-seq 表达矩阵构建与标准化

这是地基。如果这里错了，后面全崩。

1.1 定量基因表达 推荐使用 Salmon 或 Kallisto 进行拟映射（Pseudo-alignment），速度快且准确性高，最后用 tximport 将转录本水平汇总到基因水平。

# R语言示例：使用tximport汇总基因水平表达
library(tximport)
library(readr)

# 假设你有多个样本的quant.sf文件
files <- paste0("samples/", c("S1", "S2", "S3"), "/quant.sf")
names(files) <- c("S1", "S2", "S3")

# 导入转录本到基因的信息
tx2gene <- read_tsv("tx2gene.tsv") # 两列：transcript_id, gene_id

# 读取数据，countsFromAbundance="lengthScaledTPM" 是关键
# 这会将TPM按有效长度缩放，以便后续DESeq2能正确处理计数数据
txi <- tximport(files, type="salmon", tx2gene=tx2gene, 
                countsFromAbundance="lengthScaledTPM")

# 此时得到的txi$counts 可以用于DESeq2做差异分析
# txi$abundance 就是基因水平的TPM矩阵，用于后续关联分析

1.2 过滤低表达基因 不要保留所有基因。噪音太大的基因会干扰关联分析。

# 保留在至少3个样本中TPM > 1的基因
keep <- rowMeans(txi$abundance > 1) >= 0.5 # 假设50%的样本
expr_RNA <- txi$abundance[keep, ]

避坑指南 #3：处理零值（Zeros）。RNA-seq数据中，很多基因在某些样本中是0。而Log2(0)是负无穷。对于关联分析，绝对不能直接Log2。

  • 错误做法：log2(expr + 1) —— 这会导致低表达基因压缩在一起，高表达基因拉开，扭曲相关性。
  • 正确做法：使用 VST (Variable Standing Transformation) 或 rlog（在DESeq2中），或者对于质谱数据使用 MaxQuant的LFQ intensity本身（它通常已经过一定校正，但仍需Log2转换并填充缺失值）。

步骤二：蛋白质组数据预处理与缺失值填充

质谱数据最大的坑就是 Missing Values (MV)。

2.1 数据清洗 假设你从MaxQuant导出 proteinGroups.txt。

library(tidyverse)
library(imputeLCMD) # 专门处理缺失值的包

# 读取蛋白数据
prot_raw <- read.delim("proteinGroups.txt", comment.char="*")

# 提取定量列 (假设列名以 LFQ_intensity 结尾)
prot_cols <- grep("LFQ.intensity", colnames(prot_raw), value=TRUE)
prot_matrix <- prot_raw[, c("Majorityproteinidentity", prot_cols)]

# 转换为数值矩阵
prot_expr <- as.matrix(prot_matrix[, -1])
rownames(prot_expr) <- prot_matrix[, 1]
colnames(prot_expr) <- colnames(prot_matrix)[-1]

# Log2转换
prot_log2 <- log2(prot_expr)

2.2 缺失值填充（关键！） 质谱中的缺失值分为两类：

  1. MCAR (Missing Completely At Random)：随机丢失，可能是技术故障。
  2. MNAR (Missing Not At Random)：系统性的缺失，通常是因为蛋白表达量太低，低于检测限。这是最常见的情况。

如果你用均值填充，会把低表达蛋白人为抬高，产生假阳性关联。 推荐使用 left-censored imputation（左删失填充），即假设缺失值来自一个低丰度的分布。

# 使用imputeLCMD包中的QRILC方法，专门用于转录组/蛋白组数据
# 参数：k=5 (附近近邻数), seed=123 (随机种子)
prot_filled <- qrilc(prot_log2, k=5, seed=123)

# 此时prot_filled是一个填充后的矩阵，没有NA了

避坑指南 #4：不要随便用0填充或均值填充。特别是对于低丰度蛋白，均值填充会严重高估其表达量，导致与RNA数据出现虚假的高相关。务必使用专门针对 omics 数据的插补算法。

步骤三：样本匹配与批次效应校正

这一步常被忽略，但至关重要。

3.1 样本标签必须完全一致 RNA样本叫 Tumor_01，蛋白样本叫 Sample_01，这种名字对不上是找死。请确保两矩阵的列名（样本名）完全一致，并且顺序对应。

3.2 批次效应（Batch Effect） RNA和蛋白很可能是在不同时间、不同实验室、甚至不同平台测的。

  • RNA数据：用 Combat (sva包) 或 limma::removeBatchEffect 校正。
  • 蛋白数据：同样用Combat校正。

但是！关联分析前，不建议对两组数据分别做强烈的批次校正后再关联，因为校正算法可能会改变数据的分布结构，影响相关性计算。

更高级的做法是使用 ComBat-seq 或专门的整合算法（如Harmony），但在单中心研究中，通常只需确保两组数据内部已经去除了主要的批次效应。

步骤四：基因-蛋白匹配与交集获取

4.1 ID转换 RNA用的是Gene Symbol或Ensembl Gene ID，蛋白也是。确保两者使用同一套ID体系。推荐使用 Ensembl ID，因为Gene Symbol容易混淆（如旧名、同音不同义）。

4.2 取交集 只分析两组数据都检测到的基因/蛋白。

# 获取共同基因
common_genes <- intersect(rownames(prot_filled), rownames(prot_expr))

# 子集化
RNA_subset <- prot_expr[common_genes, ]
Prot_subset <- prot_filled[common_genes, ]

步骤五：关联分析与可视化

这是见证奇迹（或崩溃）的时刻。

5.1 计算相关性 对于每个基因，计算其RNA水平与蛋白水平的Pearson相关系数。

# 按行（基因）计算RNA和蛋白的相关性
cor_matrix <- cor(t(RNA_subset), t(Prot_subset), method="pearson")
gene_cor <- cor_matrix[1,] # 假设第一行是对应的蛋白，这里简化，实际需确保顺序一致

# 或者更简单的方法：直接按列对应
gene_cor <- cor(t(RNA_subset), t(Prot_subset), method="pearson")

5.2 绘制散点图（Scatter Plot） 不要只给一个相关系数，要看图。

library(ggplot2)

# 选取一个高相关的基因和一个低相关的基因做示例
gene_of_interest <- "TP53"

df <- data.frame(
  RNA = RNA_subset[gene_of_interest, ],
  Protein = Prot_subset[gene_of_interest, ],
  Condition = colnames(RNA_subset) # 如果有分组信息
)

ggplot(df, aes(x=RNA, y=Protein, color=Condition)) +
  geom_point(size=3) +
  geom_smooth(method="lm", se=FALSE, color="black") +
  theme_minimal() +
  labs(title=paste("Correlation of", gene_of_interest, ": r =", 
                   round(cor(df$RNA, df$Protein), 2)))

5.3 全局相关性分布 看看所有基因的相关性分布。

# 绘制相关性直方图
hist(gene_cor, breaks=50, main="Distribution of RNA-Protein Correlation",
     xlab="Pearson Correlation Coefficient", col="lightblue")
abline(v=0, col="red", lty=2)
abline(v=0.5, col="green", lty=2) # 标注常见阈值

第四部分：深度解读——那些“不匹配”背后的故事

当你算完相关性，一定会发现一些RNA很高但蛋白很低，或者反之的情况。这时候，新手容易停在这里，但专家会开始挖掘原因。

6.1 翻译效率分析（Translation Efficiency, TE）

TE = 蛋白丰度 / RNA丰度。

  • TE高的基因：可能受到强烈的翻译激活，或者蛋白降解慢。
  • TE低的基因：可能受到microRNA抑制，或者mRNA被隔离在P-body中。

你可以计算每个基因的TE，然后做聚类分析，看看哪些基因属于“高转录-低翻译”模块，哪些属于“转录-翻译偶联”模块。

6.2 使用Ribo-seq数据（核糖体印迹测序）

如果条件允许，加入Ribo-seq数据是解决RNA-蛋白解耦的最强工具。 Ribo-seq能直接反映正在被翻译的mRNA数量。

  • RNA-seq：总mRNA库。
  • Ribo-seq：正在被核糖体阅读的mRNA库（Polysome-associated）。
  • 蛋白组：最终产物。

如果 Ribo-seq 和 蛋白组 相关性更高，说明问题出在翻译后修饰或蛋白降解；如果 Ribo-seq 和 RNA-seq 高度一致但与蛋白组不一致，说明翻译后调控是主要因素。

6.3 基因本体（GO）富集分析

对“高相关”和“低相关”的基因分别做GO富集。

  • 高相关基因：通常参与基础代谢、核心细胞过程（如呼吸链复合物、核糖体蛋白）。这些过程往往需要快速、协调地调整蛋白量。
  • 低相关基因：通常参与信号转导、转录调控、免疫反应。这些蛋白需要“脉冲式”表达，对剂量敏感，因此转录后调控更为精细和重要。

代码示例：识别高/低相关基因并做富集

# 设定阈值，例如相关系数 > 0.7 为高相关
high_corr_genes <- names(gene_cor[gene_cor > 0.7])
low_corr_genes <- names(gene_cor[gene_cor < 0.2])

# 使用clusterProfiler做GO富集
library(clusterProfiler)
library(org.Hs.eg.db) # 人类

# 转换Gene Symbol为ENTREZ ID (GO分析需要)
high_entrez <- bitr(high_corr_genes, fromType="SYMBOL", 
                    toType="ENTREZID", OrgDb="org.Hs.eg.db")
low_entrez <- bitr(low_corr_genes, fromType="SYMBOL", 
                   toType="ENTREZID", OrgDb="org.Hs.eg.db")

# 富集分析
go_high <- enrichGO(gene = high_entrez$ENTREZID, 
                    OrgDb = org.Hs.eg.db, 
                    ont = "BP", 
                    pAdjustMethod = "BH", 
                    qvalueCutoff = 0.05,
                    readable = TRUE)

go_low <- enrichGO(gene = low_entrez$ENTREZID, 
                   OrgDb = org.Hs.eg.db, 
                   ont = "BP", 
                   pAdjustMethod = "BH", 
                   qvalueCutoff = 0.05,
                   readable = TRUE)

# 打印结果
dotplot(go_high, showCategory=10) + ggtitle("High RNA-Protein Correlation GO")
dotplot(go_low, showCategory=10) + ggtitle("Low RNA-Protein Correlation GO")

第五部分：常见坑点总结与避坑终极指南

让我把前面散落各处的坑，系统地罗列出来，方便你打印出来贴墙上。

坑点 现象 解决方案
样本名不匹配 分析报错，或者强行匹配了错误的样本 在分析前编写脚本，严格检查两组数据的样本名、顺序、数量是否完全一致。使用 intersect() 函数确保交集。
忽视Missing Values 用均值填充后，低表达蛋白假阳性关联增多 使用 QRILC 或 MinProb 等针对左删失数据的填充方法。明确报告缺失值比例。
直接Log2(0) 报错或产生极端离群值 先填充缺失值，再进行Log2转换。对于RNA-seq，使用VST变换优于简单的Log2(TPM+1)。
混合不同标准化方法 一组用TPM，一组用Counts，导致尺度完全不同 关联分析前，确保两组数据都在同一尺度（通常是对数刻度下的标准化值）。TPM和LFQ Intensity都可以Log2后直接比较。
忽略基因异构体 一个基因有多个蛋白亚型，质谱只打到其中一条肽段 尽量在基因水平进行关联。如果必须做异构体水平，需要确保质谱数据有异构体特异性肽段（Isoform-specific peptides），否则结果不可信。
批次效应混淆 差异表达主要来源于测序
-- 展开阅读全文 --