第4章 蛋白质组学数据基础：质谱数据、肽段鉴定、蛋白质丰度矩阵
聊到蛋白质组学，很多人第一反应是“啊，质谱”。没错，质谱就是蛋白质组学最主力的数据来源。但说实话，我刚入行时对着那堆 m/z 和强度值完全懵了——这玩意儿怎么能变成深度学习能用的张量？这一章咱们就把这层窗户纸捅破。
4.1 质谱数据长什么样？
质谱的本质很简单：把肽离子按质荷比（m/z）分开，记录每个 m/z 处有多少离子（强度）。输出就是一个二维谱图，横轴 m/z，纵轴强度。
实际项目中你遇到的往往是 RAW 或者 mzML 格式。我建议初学者先别碰二进制格式，直接从 mzXML 或 mzML 的文本化版本入手。下面是常见格式对比：
| 格式 | 特点 | 典型处理工具 |
|---|---|---|
| RAW（Thermo） | 二进制，厂商锁定 | Thermo RawFileReader |
| mzML | HUPO-PSI 标准，XML | pymzML, ProteoWizard |
| mzXML | 早期开放格式 | Trans-Proteomic Pipeline |
| MGF | 纯文本谱图列表 | 几乎任何工具 |
我记得第一次做大规模分析时，数据量 50GB 起。你不能把整个质谱原始数据全塞进 GPU——那几百亿个数据点多半是噪音。咱们要的是经过提取的肽段谱图（MS2）和对应的定量信息。
小技巧：用 ProteoWizard 的 msConvert 可以把 RAW 转成 mzML，同时做 peak picking。命令行如下（别复制错了）：
msconvert --mzML --filter "peakPicking true 1-" input.raw -o output_dir加 --zlib 可以压缩大小，我一般会加。
4.2 肽段鉴定——从谱图到序列
质谱仪打到肽段上，产生碎片离子。肽段鉴定就是拿这些碎片模式去和理论谱图匹配。最常见的算法是 SEQUEST、Mascot、Andromeda（MaxQuant 用的），还有现在的深度学习工具比如 DeepMass、Prosit。
为什么我单独提这段？因为很多组学模型训练的数据集就是直接拿鉴定结果当标签。如果鉴定步骤本身有错，模型学到的就是天书。
我曾在一个项目里遇到过这件事：某公共数据库的肽段鉴定 FDR 设到了 5%，结果模型训练完看起来 AUPRC 很高，但实际预测全是假阳性。搞得我排查了两周。从那以后我经手的项目，鉴定 FDR 必须 ≤ 1%，最好 0.5%。
鉴定结果通常存成 pepXML 或 mzid 格式，核心字段就是：
- 扫描号（Scan）
- 肽段序列
- 电荷状态
- 置信度分数（比如 XCorr, e-value, q-value）
4.3 蛋白质丰度矩阵——你的训练目标
肽段鉴定后，接下来是把肽段映射到蛋白质，并算出每个蛋白质在各样本里的丰度。这就得到了丰度矩阵：行是蛋白质，列是样本，单元格里是丰度值（通常是对数转换后的 TMT 报告离子强度或 LFQ intensity）。
这个矩阵就是深度学习模型的核心训练数据了。
- 输入：蛋白质序列特征（比如 embedding）或其他组学特征
- 输出：这个蛋白质在疾病/正常条件下的差异丰度
- 注意：矩阵里缺失值非常多——质谱有天然“缺失采样”问题
避坑指南：不要直接把原始强度值丢给 GPU。我试过不归一化就训练，Loss 下不去。必做步骤：
- log2 转换（让数据接近正态）
- 分位数归一化（消除批次效应）
- 缺失值处理（我倾向于用 KNN imputation 而非填零）
4.4 数据组织与输入管道
拿到丰度矩阵后，为了喂给 GPU，你得设计一个能高效读取的数据管道。这里我分享一个在实践中打磨过的策略：
# 以 PyTorch 为例，但思路通用
dataset = ProteomicsDataset(protein_ids, expression_data, sample_labels)
dataloader = DataLoader(dataset, batch_size=32,
num_workers=8, pin_memory=True)
但千万别把整个矩阵都加载进内存再 shuffle——尤其当你处理的是 10 万蛋白 × 1000 样本时，内存直接爆掉。我一般会做一个内存映射文件（memory-mapped）或者用 h5py 按批次读取。
# 使用 h5py 存储丰度矩阵
import h5py
f = h5py.File('proteomics.h5', 'w')
f.create_dataset('expression', shape=(n_proteins, n_samples),
dtype='float32', compression='gzip')
f.create_dataset('protein_ids', shape=(n_proteins,), dtype='S20')
实际加载时只 select 当前 batch 需要的行，配合 GPU 的 pinned memory 做异步传输——这样可以掩盖 IO 延迟。
你自己想想：一次推理一个样本和一次推理一个 batch，训练速度差多少？我见过有人单样本循环跑了一天一夜还没跑完，改成 batch 后两小时搞定。
4.5 全流程可视化：从质谱原始输出到丰度矩阵
下面我用一张图把整个数据流串起来。这也是我每次给新同学讲时必定先画的框架。
4.6 实际项目中的常见坑
- 蛋白质 ID 不一致：不同数据库（UniProt、RefSeq）ID 格式不同，合并前必须统一映射。我有个脚本专门做这件事。
- 肽段到蛋白质的归纳（protein inference）：同源肽段会被多个蛋白质共享，这叫“ambiguously assigned peptide”。我建议删除这些共享肽段，不然会引入噪声。
- 批次效应：TMT 多批次实验之间的丰度漂移非常明显。用 ComBat 或 limma 调整是必须的。
核心总结：质谱数据经过谱图预处理、肽段鉴定、蛋白推断三步，才能得到可训练的丰度矩阵。矩阵里缺失多、噪声大，但 GPU 训练的前提是先把数据组织好——管道设计比模型设计更容易被忽略。把这一步做扎实，模型训练就成功了一半。