质谱蛋白组学数据分析¶
一句话概述¶
利用MaxQuant进行DDA数据蛋白鉴定与LFQ定量、DIA-NN处理DIA数据实现深度覆盖蛋白组定量，再用Perseus进行统计分析和生物学解读，是蛋白质组学研究的标准计算流程。
核心知识点表格¶
| 知识点 | 关键内容 | 常用工具 |
|---|---|---|
| 数据采集模式 | DDA(数据依赖)vs DIA(数据独立) | - |
| 数据库搜索 | 肽段-谱图匹配(PSM) | MaxQuant/Andromeda, MSFragger |
| FDR控制 | Target-Decoy策略，1% FDR | Percolator, PeptideProphet |
| 定量方法 | LFQ/TMT/iTRAQ/SILAC | MaxQuant, Proteome Discoverer |
| DIA分析 | 谱库搜索或library-free | DIA-NN, Spectronaut, OpenSWATH |
| 统计分析 | 缺失值填充、差异检验、聚类 | Perseus, limma, MSstats |
| 功能注释 | GO/KEGG/蛋白互作网络 | STRING, g:Profiler |
各步骤详解¶
第一步：理解质谱蛋白组学基本原理¶
白话解释： 蛋白质组学的基本流程（Bottom-up）：蛋白质→酶解为肽段→液相色谱分离→质谱检测→计算机搜库鉴定。质谱仪测量的是肽段的质荷比(m/z)和碎片谱图，通过与蛋白质数据库中的理论谱图比较来鉴定肽段身份。
技术细节： - DDA (Data Dependent Acquisition): 选择Top-N强度离子进行MS2碎裂 - DIA (Data Independent Acquisition): 对所有m/z窗口系统性碎裂 - DDA优势：鉴定可靠、搜库简单；劣势：随机采样、重复性差 - DIA优势：覆盖度高、重复性好；劣势：数据解析复杂 - 典型深度：DDA ~5000-8000蛋白，DIA ~8000-12000蛋白
第二步：MaxQuant DDA数据分析¶
白话解释： MaxQuant是DDA数据分析的金标准工具。它把原始质谱文件（.raw）与蛋白质数据库比对，找出每个谱图对应哪个肽段（PSM），然后把肽段归属到蛋白质，最后用LFQ算法进行无标记定量。
技术细节： - Andromeda搜索引擎：基于概率打分的PSM匹配 - MaxLFQ算法：基于肽段比值的蛋白定量，不受总量影响 - Match Between Runs (MBR)：利用保留时间对齐补全缺失值 - Razor peptide：共享肽段归属给最可能的蛋白质组
代码示例：
# MaxQuant通常通过GUI运行
# 但也支持命令行模式
# 准备参数文件 mqpar.xml（从GUI导出或手动编写）
# 关键参数设置：
# 命令行运行MaxQuant
dotnet MaxQuant/bin/MaxQuantCmd.exe mqpar.xml
# 或Linux下：
# MaxQuant 2.x 已切换为 .NET Core，Linux 下可直接用 dotnet MaxQuantCmd.dll
mono MaxQuant/bin/MaxQuantCmd.exe mqpar.xml
MaxQuant关键参数配置：
<!-- mqpar.xml 核心参数 -->
<MaxQuantParams>
<!-- 数据库 -->
<fastaFiles>
<string>uniprot_human_reviewed_2024.fasta</string>
<string>contaminants.fasta</string>
</fastaFiles>
<!-- 酶切 -->
<enzymes>
<string>Trypsin/P</string>
</enzymes>
<maxMissedCleavages>2</maxMissedCleavages>
<!-- 修饰 -->
<fixedModifications>
<string>Carbamidomethyl (C)</string>
</fixedModifications>
<variableModifications>
<string>Oxidation (M)</string>
<string>Acetyl (Protein N-term)</string>
</variableModifications>
<!-- 定量 -->
<lfqMode>1</lfqMode> <!-- LFQ开启 -->
<lfqMinRatioCount>2</lfqMinRatioCount>
<matchBetweenRuns>True</matchBetweenRuns>
<matchTimeWindow>0.7</matchTimeWindow>
<!-- FDR -->
<proteinFdr>0.01</proteinFdr>
<peptideFdr>0.01</peptideFdr>
<siteFdr>0.01</siteFdr>
</MaxQuantParams>
# MaxQuant输出文件（在combined/txt/目录下）：
# proteinGroups.txt - 蛋白质定量结果（主要分析文件）
# peptides.txt - 肽段级别结果
# evidence.txt - PSM级别详细信息
# msms.txt - 所有MS2谱图
# summary.txt - 运行摘要统计
# parameters.txt - 使用的参数
# 质量检查关键指标：
# - MS/MS identification rate: >25% (理想>40%)
# - Peptide FDR: ≤1%
# - Protein FDR: ≤1%
# - 蛋白质鉴定数: human cell line通常>5000
第三步：DIA-NN数据分析¶
白话解释： DIA-NN是处理DIA数据最快且性能优秀的工具之一。它可以使用预建谱库（library-based）或直接从DIA数据中"自学"鉴定肽段（library-free模式）。后者不需要额外的DDA数据建库，更加便捷。
技术细节： - Library-free模式：从FASTA直接预测肽段碎片谱+保留时间 - 深度学习预测：用神经网络预测MS2谱图和RT - Double-pass搜索：第一轮建立经验性谱库，第二轮精确定量 - MaxLFQ归一化集成 - 典型速度：~10分钟/raw文件
代码示例：
# 安装DIA-NN
# 下载: https://github.com/vdemichev/DiaNN
# ⚠ DIA-NN 2.0（2024）重大变更：library-free 模式需两步运行
# 步骤1：生成 in silico 谱库（不加载 raw 文件）
# 步骤2：加载 raw 文件和步骤1生成的谱库进行分析
# 以下命令适用于 DIA-NN 1.x；使用 2.0 请参考官方文档
# Library-free模式（推荐，无需DDA谱库）
diann --f raw_files/ \
--lib "" \ # 空=library-free模式
--fasta uniprot_human.fasta \
--fasta-search \
--out diann_output/report.tsv \
--qvalue 0.01 \
--matrices \ # 输出定量矩阵
--threads 16 \
--predictor \ # 使用深度学习预测
--gen-spec-lib \ # 生成实验特异性谱库
--smart-profiling \
--met-excision \ # N-端Met切除
--cut K*,R* \ # Trypsin切割规则
--missed-cleavages 2 \
--min-pep-len 7 \
--max-pep-len 30 \
--min-pr-charge 2 \
--max-pr-charge 4 \
--var-mods 1 \ # 可变修饰数
--var-mod UniMod:35,15.994915,M \ # M氧化
--unimod4 # C碘乙酰胺化
# Library-based模式（使用预建谱库）
diann --f raw_files/ \
--lib spectral_library.tsv \
--out diann_output/report.tsv \
--qvalue 0.01 \
--matrices \
--threads 16
# 输出文件：
# DIA-NN 2.0+ 输出格式从 .tsv 变更为 .parquet，使用 pandas/pyarrow 读取
# report.tsv - 主要结果（肽段/蛋白鉴定+定量）
# report.pr_matrix.tsv - 前体离子定量矩阵
# report.pg_matrix.tsv - 蛋白质组定量矩阵（MaxLFQ）
# report.stats.tsv - 每个run的统计信息
第四步：Perseus统计分析¶
白话解释： Perseus是MaxQuant配套的下游统计分析平台。主要工作：1）过滤掉污染物和反库蛋白；2）处理缺失值；3）做差异表达蛋白检测（t-test或ANOVA）；4）聚类和可视化。
技术细节： - 缺失值类型：MNAR(缺失不随机，低丰度蛋白检测不到) vs MAR(随机缺失) - MNAR用左移正态分布填充（模拟低于检测限的值） - 有效过滤：至少在一组中70%的样本有值 - 多重检验校正：permutation-based FDR (s0参数) - s0参数：控制fold change和p-value的权衡
代码示例：
# 使用R代替Perseus进行统计分析（更灵活）
library(tidyverse)
library(limma)
library(impute)
library(pheatmap)
library(EnhancedVolcano)
# 1. 读取MaxQuant/DIA-NN输出
pg <- read.table("proteinGroups.txt", sep="\t", header=TRUE,
quote="", comment.char="")
# 2. 过滤
pg_clean <- pg %>%
filter(Reverse != "+") %>% # 去除反库蛋白
filter(Potential.contaminant != "+") %>% # 去除污染物
filter(Only.identified.by.site != "+") # 只靠修饰位点鉴定的去除
# 3. 提取LFQ intensity矩阵
lfq_cols <- grep("^LFQ.intensity.", names(pg_clean), value=TRUE)
lfq_matrix <- pg_clean[, lfq_cols]
rownames(lfq_matrix) <- pg_clean$Majority.protein.IDs
# 4. Log2转化
lfq_log2 <- log2(lfq_matrix)
lfq_log2[lfq_log2 == -Inf] <- NA
# 5. 有效值过滤（至少在一组中有70%的有效值）
groups <- c("ctrl","ctrl","ctrl","treat","treat","treat")
valid_in_group <- function(x, group_labels) {
any(tapply(!is.na(x), group_labels, function(g) mean(g) >= 0.7))
}
keep <- apply(lfq_log2, 1, valid_in_group, group_labels=groups)
lfq_filtered <- lfq_log2[keep, ]
# 6. 缺失值填充（MNAR: 左移正态分布）
impute_mnar <- function(mat, shift=1.8, width=0.3) {
mat_imp <- mat
for(i in 1:ncol(mat)) {
valid_values <- mat[!is.na(mat[,i]), i]
mu <- mean(valid_values) - shift * sd(valid_values)
sigma <- width * sd(valid_values)
n_missing <- sum(is.na(mat[,i]))
mat_imp[is.na(mat[,i]), i] <- rnorm(n_missing, mu, sigma)
}
mat_imp
}
lfq_imputed <- impute_mnar(as.matrix(lfq_filtered))
# 7. 差异表达分析（limma）
design <- model.matrix(~0 + factor(groups))
colnames(design) <- c("ctrl", "treat")
contrast <- makeContrasts(treat - ctrl, levels=design)
fit <- lmFit(lfq_imputed, design)
fit2 <- contrasts.fit(fit, contrast)
fit3 <- eBayes(fit2)
results <- topTable(fit3, number=Inf, adjust.method="BH")
# 8. 筛选差异蛋白
sig_proteins <- results %>%
filter(adj.P.Val < 0.05, abs(logFC) > 1)
cat("上调蛋白:", sum(sig_proteins$logFC > 0), "\n")
cat("下调蛋白:", sum(sig_proteins$logFC < 0), "\n")
# 9. 火山图
EnhancedVolcano(results,
lab = rownames(results),
x = 'logFC', y = 'adj.P.Val',
pCutoff = 0.05, FCcutoff = 1,
title = "Differential Protein Expression")
第五步：功能富集与蛋白互作网络¶
白话解释： 找到差异蛋白后，用GO/KEGG富集看这些蛋白参与什么生物过程，用STRING数据库看它们之间有没有已知的蛋白-蛋白互作关系。
代码示例：
# GO/KEGG富集
library(clusterProfiler)
library(org.Hs.eg.db)
# 转换蛋白ID为基因名
gene_list <- gsub(".*\\|(.*)\\|.*", "\\1", rownames(sig_proteins))
ego <- enrichGO(gene = gene_list,
OrgDb = org.Hs.eg.db,
keyType = "UNIPROT",
ont = "BP",
pAdjustMethod = "BH")
# STRING蛋白互作网络
# 使用STRINGdb R包
library(STRINGdb)
string_db <- STRINGdb$new(version="12.0", species=9606, score_threshold=400)
mapped <- string_db$map(sig_proteins_df, "gene_name", removeUnmappedRows=TRUE)
string_db$plot_network(mapped$STRING_id[1:50])
# 导出用Cytoscape可视化
string_db$get_interactions(mapped$STRING_id) %>%
write.table("string_interactions.tsv", sep="\t", row.names=FALSE)
第六步：MSstats标准化定量分析¶
白话解释： MSstats是专门为蛋白质组学设计的统计框架，能更好地处理蛋白质组学数据的特殊性质（如缺失值模式、异方差性），支持复杂实验设计。
代码示例：
library(MSstats)
# 从MaxQuant结果转换
# 需要evidence.txt和proteinGroups.txt
raw <- MaxQtoMSstatsFormat(
evidence = read.table("evidence.txt", sep="\t", header=TRUE),
proteinGroups = read.table("proteinGroups.txt", sep="\t", header=TRUE),
annotation = read.csv("annotation.csv"), # Run, Condition, BioReplicate
removeProtein_with1Feature = TRUE)
# 数据处理（归一化+汇总）
processed <- dataProcess(raw,
normalization = "equalizeMedians",
summaryMethod = "TMP", # Tukey's Median Polish
censoredInt = "NA",
MBimpute = TRUE, # 基于模型的缺失值填充
maxQuantileforCensored = 0.999)
# 差异分析
comparison <- matrix(c(-1, 1), nrow=1)
colnames(comparison) <- c("Control", "Treatment")
rownames(comparison) <- "Treatment-Control"
test_result <- groupComparison(
contrast.matrix = comparison,
data = processed)
# 结果
sig_msstats <- test_result$ComparisonResult %>%
filter(adj.pvalue < 0.05, abs(log2FC) > 1)
# 可视化
groupComparisonPlots(test_result$ComparisonResult,
type = "VolcanoPlot",
sig = 0.05, FCcutoff = 1)
实战命令¶
#!/bin/bash
# === 蛋白质组学分析完整流程 ===
# DIA-NN library-free分析
diann --f /data/raw/*.raw \
--lib "" --fasta human_uniprot.fasta --fasta-search \
--out results/report.tsv --qvalue 0.01 --matrices \
--threads 32 --predictor --gen-spec-lib --smart-profiling \
--met-excision --cut K*,R* --missed-cleavages 2 \
--var-mods 1 --var-mod UniMod:35,15.994915,M --unimod4
# R下游分析
Rscript proteomics_analysis.R
面试常问点¶
Q1: DDA和DIA的区别是什么？ A: DDA选择当前扫描中最强的N个离子进行MS2碎裂（数据依赖），导致低丰度蛋白随机遗漏；DIA对预设的m/z窗口系统性碎裂所有离子，理论上不遗漏，但谱图更复杂需要专门算法解析。DIA重复性更好、覆盖更深，是当前趋势。
Q2: FDR 1%是如何计算的？ A: Target-Decoy策略：将蛋白序列反转(reversed)构建诱饵数据库，与目标数据库同时搜索。假设诱饵库的匹配全是假阳性，FDR = 2×诱饵匹配数/总匹配数（在某打分阈值下）。蛋白水平FDR通过parsimony原则从肽段FDR推断。
Q3: LFQ定量的原理？ A: MaxLFQ算法基于肽段pair-wise比值的中位数估计蛋白质相对丰度。不需要标记，直接利用MS1信号强度。优势：简单灵活；劣势：缺失值多，需要较多样本。
Q4: 蛋白组学缺失值如何处理？ A: 区分两类缺失：MNAR（低于检测限，用左移正态分布填充）和MAR（随机缺失，用KNN/SVD等方法填充）。实践中常用策略：先过滤（至少在一组中70%样本有值），再用混合方法填充。也可使用基于模型的方法（如MSstats的TMP+Imputation）。
Q5: TMT和LFQ的选择？ A: TMT：同位素标记，最多18plex，技术变异小，但ratio compression问题限制动态范围，成本高。LFQ：无标记，灵活便宜，但技术变异大需要更多重复。大规模队列研究推荐DIA-LFQ；精确定量少量样本推荐TMT。
易错点¶
数据库版本不匹配：使用过时或错误物种的数据库会严重影响鉴定率。务必使用最新UniProt reviewed数据库+污染物库。
缺失值处理不当：全部用零替代或不区分MNAR/MAR类型的填充会引入系统性偏差。
蛋白推断歧义：共享肽段导致蛋白鉴定不确定性。报告时使用protein groups而非单个accession。
忽略批次效应：多批次TMT实验需要桥接通道(bridge channel)校正批次效应。
Match Between Runs过度填充：MBR可能在低质量样本中填充不可靠的值，需要评估其质量。
统计方法选择不当：蛋白组学数据不满足标准t-test假设，推荐使用limma(经验贝叶斯)或MSstats等专用方法。
补充知识¶
蛋白质组学搜索引擎比较¶
| 搜索引擎 | 速度 | 特点 | 适用场景 |
|---|---|---|---|
| MaxQuant/Andromeda | 慢 | LFQ金标准、MBR | DDA-LFQ |
| MSFragger | 极快 | open search、glyco | 大数据集/修饰发现 |
| DIA-NN | 快 | 深度学习预测谱库 | DIA |
| Spectronaut | 中 | 商业，hybrid library | DIA(高端) |
| Proteome Discoverer | 中 | 商业，Thermo整合 | TMT/DDA |