xcms这个软件说实话已经快20岁了。一个在2005年推出的R包到现在依旧是代谢组学数据处理绕不开的选项这本身就挺罕见的。我早期做非靶向代谢组学的时候就被它折腾过也被它惊艳过。今天想系统的聊聊这个“老伙计”聊聊它到底解决了什么问题以及为什么到了今天我们处理LC-MS数据时还是要把它的逻辑摸清楚。1. 核心设计与定位为什么2005年需要xcms现在依然需要1.1 当年的困境数据处理的“黑暗森林”在xcms出现之前做代谢组学或者相关质谱数据分析是个非常痛苦的事情。那个年代仪器厂商各自为政输出格式五花八门数据分析基本靠仪器自带的工作站配合一堆手工Excel操作。你想把一个样本的原始数据批量处理或者想换个算法重新处理一遍几乎等于噩梦。更麻烦的是那时候高分辨质谱开始普及数据量猛增。一台Q-TOF或者Orbitrap跑完一批样本产生的原始文件动辄几个G。传统的手动积分、逐峰确认方式根本没法应付。代谢组学的核心是“非靶向”意味着你不知道你要找什么你只知道一堆样本里有差异。这个“不知道”就需要一个能自动化、能重现、能处理大规模数据的工作流。xcms的价值就是在这个时间点把整个LC-MS数据处理流程——从原始文件读入、峰检测、峰对齐到最终输出一个表达矩阵——用R语言这套开放生态变成了一条标准化的流水线。它借了Bioconductor的东风让统计分析PCA、OPLS-DA等和数据处理在同一个环境里无缝衔接这对于当时习惯了“工作站处理完导Excel再进SPSS”的人来说是降维打击。1.2 解题思路“峰表”为核心替代手工积分xcms的整体设计思路用一句话概括就是先把复杂的三维原始数据保留时间、质荷比、强度抽象成一张“峰表”之后的统计分析全都基于这张表展开。这个思路现在看着顺理成章但在当时是很有魄力的。因为质谱原始数据是连续的谱图而代谢组学关心的是“哪些化合物在哪些样本里有多少”。xcms通过算法把信噪比符合要求的信号峰拉出来记录它们的m/z、保留时间、峰面积或峰高然后通过保留时间对齐和峰匹配把不同样本里同一个化合物对应起来最后输出一张“样本 × 峰”的表格。这个设计的优势太明显了它在数据处理和分析之间画了一条清晰的界线。数据处理解决“这个峰是什么、有多大”统计分析解决“这些峰有没有统计学意义”。你可以在任意阶段停下来检查也可以换不同的算法参数重新跑而整个流程是可重复、可记录的。这一点直到今天很多商业软件都做得不如xcms彻底。2. 核心算法拆解centWave到底强在哪里2.1 从matchedFilter到centWave只要是老用户都知道xcms最初的默认峰检测算法其实不是centWave而是matchedFilter。这个算法的逻辑有点像“模板匹配”它在一系列的m/z窗口里对信号做平滑处理然后设定一个阈值来找“小山包”。但matchedFilter有个致命弱点它对峰宽的容忍度很低。如果你的色谱峰比较宽或者基线漂移严重它就容易漏检或者把噪声当成峰。而且参数不好调一个数据集调好的参数换一批样本又不行了。好在xcms在2006年左右引入了centWave算法这才是真正让它“正中代谢组学眉心”的大杀器。centWave的思路完全不同它不依赖固定的m/z窗口模板而是先在原始数据里寻找“感兴趣区域”ROIRegions of Interest。它的逻辑是一个真正的色谱峰在连续的质谱扫描中它的精确m/z应该是基本不变的变化的只是强度。所以算法先在m/z维度上扫描把所有能连续匹配上且强度变化符合色谱峰形状的信号串联起来形成ROI候选。然后再对这些候选的色谱峰做连续小波变换CWTContinuous Wavelet Transform通过小波变换的多尺度特性识别峰的真实形状和边界。2.2 用生活的类比理解小波变换这个连续小波变换听起来抽象可以这样理解色谱峰的形状就像一座起伏的山丘但山丘大小不一高的矮的宽的窄的都有。你用一个固定大小的尺子去量肯定量不准。小波变换相当于拿了一把可以自动伸缩的尺子——它从多个尺度去扫描信号哪个尺度下信号和小波的形状最匹配就认为那里有个峰并顺手把峰的宽度和位置定下来。这个设计的聪明之处在于它对重叠峰和低丰度信号的识别能力。实际样本里杂质非常多两个峰挨得很近是常态。centWave先通过ROI锁定候选再用小波变换精确定位这套组合拳在抗干扰性上远比简单的滑动窗口平滑要强。我现在用的很多新软件算法内核还是在吃centWave的红利。3. 实操工作流从原始数据到可统计的峰表3.1 数据准备格式转换别给自己埋雷无论用哪个工具第一步永远是数据格式转换。不要用仪器厂商自己的格式直接喂给xcms除非你想体验掉头发的感觉。厂商格式.wiff、.raw、.d等往往封装了私有信息甚至有些还带着加密xcms读起来又慢又容易出错。我的习惯是先把所有原始文件用msConvert (ProteoWizard的一部分) 统一转成.mzML格式。转换时有几个细节要注意profile data和centroid data的选择早期很多教程建议用centroid数据质心模式因为文件小、处理快。但如果你用的是centWave我强烈建议用profile数据轮廓模式来跑因为centWave在识别ROI时会用到谱峰的形状信息profile模式结果更稳。虽然文件大一点但是为了准确性值。文件命名要规范尽量短且不含特殊字符不要用中文、不要带空格更不要用“-”号。R在读取文件名时遇到这些容易出问题而且后续写代码做分组时文件名要用来批量提取样本组信息。推荐格式sample1_control_rep1.mzML、sample2_disease_rep2.mzML这种用下划线分好信息段后面一个strsplit就能拿到分组信息。保留时间、m/z校准转格式时不要乱动数据关闭所有降噪、平滑选项保持“原汁原味”。xcms的算法是设计来处理原始噪声的你提前平滑反而可能引入假峰或者破坏低丰度信号的真实形状。3.2 读入数据封装成XcmsExperiment对象现在的xcms是第三版架构了读入数据的方式和老的xcmsSet对象相比变化巨大。现在的统一入口是readMSData或者直接用xcms::xcmsExperiment函数。读入后数据会被封装成一个XcmsExperiment对象后续所有操作都围绕这个对象转。# 建议在R 4.2以上版本运行安装最新版xcms library(xcms) library(MSnbase) # 读取所有mzML文件注意用phenoData描述样本分组信息 raw_data - readMSData( files list.files(your_path, pattern .mzML, full.names TRUE), mode onDisk ) # 查看基本信息 raw_data如果你的样本分成了两组一定要在读取之后给每个样本补充分组信息。这一步很多人会忽略导致后面做PCA或者差异分析时不知道怎么分组。正确的做法是这样# 假设文件名都是 sample_group_rep.mzML 格式 sample_meta - data.frame( sample_name basename(fileNames(raw_data)), stringsAsFactors FALSE ) sample_meta$group - gsub(.*_([a-zA-Z])_.*, \\1, sample_meta$sample_name) # 把分组信息塞进对象 sampleData(raw_data)$group - sample_meta$group3.3 峰检测参数调优的艺术核心中的核心峰检测是整个流程中最关键的步骤参数选错了后面对齐、注释做得再漂亮也是白搭。我用的是findChromPeaks函数算法选CentWaveParam。# 定义centWave参数 cwp - CentWaveParam( peakwidth c(5, 20), # 色谱峰宽范围秒 snthresh 10, # 信噪比阈值 noise 1000, # 噪声过滤阈值 ppm 15, # 质量精度偏差ppm mzcenter TRUE, integrate 1, # 峰积分方法 mzdiff -0.001, # 相邻峰最小m/z差 fitgauss TRUE, # 尝试高斯拟合 verboseColumns TRUE # 输出额外诊断信息 ) # 执行峰检测 xdata - findChromPeaks(raw_data, param cwp)这些参数里最需要花心思的是peakwidth和ppm。peakwidth必须根据你的色谱条件设置。如果你用的是UPLC峰宽可能集中在5-15秒如果是普通HPLC峰宽可能到20-40秒。不要偷懒用默认值最合理的办法是打开一个代表性样本的TIC图随机挑几个明显的峰量一下它们的峰底宽和峰顶宽再填进去。ppm则取决于你的质谱仪精度。Orbitrap可能只需要5-10Q-TOF可能给15-25给得太小会漏峰给得太大容易把同位素峰也当成候选增加后续去重负担。snthresh设为10就够用了。有时候为了尽量召回真峰可以先设到5后面靠统计手段过滤。但如果你直接设到30或者更高很可能把低丰度的代谢物全丢了做非靶向这是不能接受的。我一般先设810跑一轮看峰数量曲线再做调整。记住峰检测宁可多召回不可漏检。后续可以用peakQC等工具或者用QC样本的变异系数CV来过滤不可靠峰。3.4 分组对齐跨样本匹配决定数据质量检测完峰每个样本里都有一堆峰表。但不同样本之间同一个化合物会因为梯度漂移、温度变化等原因保留时间发生偏移。你需要让算法知道“样本1里RT120.5s的峰”和“样本2里RT123.1s的峰”其实是同一个东西。这就是对齐分组。# 先用peak density方法做初始分组 pdp - PeakDensityParam( sampleGroups sampleData(xdata)$group, bw 5, # 带宽允许的保留时间漂移范围秒 minFraction 0.5, # 该峰在至少50%的样本中存在 binSize 0.025 # m/z方向的bin宽度 ) xdata - groupChromPeaks(xdata, param pdp) # 再进行保留时间校正 xdata - adjustRtime(xdata, param ObiwarpParam())bw参数是核心它控制聚类时允许的保留时间漂移范围。如果你的样本批次比较稳定bw2或3就够了如果样本量大或者跨度时间长、漂移严重设置5甚至10都不过分。但要注意bw设得太大会把本来不同的峰错误地匹配在一起。最关键的是这一步做完一定要检查。绘制保留时间校正前后的偏差图看看校正曲线是否平滑有没有明显的跳动。如果校正后反而更乱了就需要考虑是不是binSize设得太大或者太小。对齐之后再做一次分组第二次分组这一次就是最终确定峰表了。有些教程会省略第二次分组直接fillChromPeaks以前我这么干过后来发现峰表里有很多NA就是因为最终分组没有做。正确操作是# 重新分一次组生成最终的峰表 xdata - groupChromPeaks(xdata, param pdp) # 填充缺失峰用积分区域内的原始信号估算缺失值 xdata - fillChromPeaks(xdata, param FillChromPeaksParam())fillChromPeaks这一步很关键。一张表达矩阵里全是NA下游统计做着很费劲。这个函数会用已知峰的保留时间和m/z窗口去原始数据里重新积分把缺失值补上。注意它不等于“造数据”补进去的值是有原始信号支撑的这样后续做PCA和差异分析时会顺利很多。3.5 提取结果从峰表到数据矩阵上面几步跑完就可以从xdata对象里提取最核心的结果了——那个用于后续分析的峰表。# 获取峰表 peaks_tbl - chromPeaks(xdata)如果你想把峰面积变成“行是峰列是样本”的矩阵直接用featureValues# 提取特征表达矩阵 feature_mat - featureValues( xdata, method maxint, # 同个特征在样本中出现多个峰时取强度最大者 value into, # 使用峰面积 filled TRUE # 包含填充值 ) # 转换成数据框并导出 feature_df - as.data.frame(feature_mat) write.csv(feature_df, xcms_feature_matrix.csv, quote FALSE)这里的value参数值得多说两句。into是峰面积maxo是峰高。大部分代谢组学文献用峰面积作半定量如果你的色谱峰形普遍不好比如严重拖尾用峰高反而可能更稳定。我一般两个都会试比较一下QC样本CV值再决定用哪个出结果。4. 进阶技巧跑批正确别只盯着默认值4.1 用IPO自动优化参数但别迷信自动化xcms最劝退新手的地方就是参数太多调起来全靠手感。有个R包叫IPOIsotopic Peak Optimization能做参数自动调优。它的原理是跑多组参数组合用一个响应函数评估峰数量、峰形、重复性等指标的平衡。我用下来觉得IPO作为初筛工具很好用尤其是你拿到一个新类型的数据、完全没概念时跑一遍IPO能拿到一套不错的起步参数。library(IPO) # 这是早期版本用法新版本请参考IPO文档 peaks_ipo - findPeaks_IPO( raw_data, params list(ppm c(10, 20), snthresh c(5, 10), peakwidth c(5, 20)) )但IPO也有局限。第一它很费时间大数据集跑一个组合要半小时全部跑完可能一天就没了。第二它优化的是“统计意义上的最优”未必是“生物学意义上的最优”。它可能会为了减少峰数量而牺牲掉一批低丰度但很有意义的小峰。所以我自己习惯的做法是先用IPO拿初始参数然后人工抽查几张总离子流图看看那些你关心的目标峰、或者已知的内标峰有没有被正确检出再微调noise和snthresh。4.2 QC样本和空白的处理方式拿到峰表之后QC样本就是你的照妖镜。做非靶向代谢组学的标准流程隔几个样本就穿插一个QC样本同一个混合样本分装成很多针。这批QC样本要参与先前的峰检测和对齐但在最终统计中它们用来评估技术重复性。计算每个特征的QC样本相对标准偏差RSD凡是RSD 30%的峰原则上都应该剔除因为这样的峰不可靠。# 假设QC样本名里有QC字样 qc_idx - which(grepl(QC, colnames(feature_mat))) rsd - apply(feature_mat, 1, function(x) { sd(x[qc_idx], na.rm TRUE) / mean(x[qc_idx], na.rm TRUE) }) keep - which(rsd 0.3) filtered_mat - feature_mat[keep, ]空白样本blank也很有用。做样本前先跑一针空白溶剂这一步不是为了找峰而是为了记录背景信号。后面过滤特征时凡是在空白里也有明显信号的峰大概率是溶剂残留、塑料增塑剂或者仪器污染要剔除。4.3 把峰表转换成能做差异分析的格式xcms输出的峰表只是一个宽表“行是特征列是样本”。要真正做差异代谢物筛选通常还要做缺失值填补、数据归一化、Pareto缩放这些步骤然后再扔给ropls做OPLS-DA。我习惯在xcms峰表基础上用metabolomics这类工具做后续处理。但无论如何xcms这一步输出的峰表质量直接决定了后面能筛出来多少靠谱的差异代谢物。很多人在这一步草草了事后面拼命调统计参数其实是本末倒置了。5. 常见问题与排查技巧你大概率会踩的这些坑5.1 峰数量少得可怜或者多到爆炸跑完findChromPeaks发现检测出的峰数量和你预期严重不符。这时候不要慌按顺序排查峰太少第一反应是snthresh设太高了把它从10降到5试试。第二个可能是ppm设太小导致ROI识别时很多谱峰没有被连接起来。第三个可能是noise假设如果质谱仪器基线噪声很高设一个高于实际噪声水平的noise值会直接把低丰度峰全屏蔽掉了。峰太多多半是peakwidth范围给太宽了很多噪声也被当成了峰。比如你把peakwidth下限设成1秒那几乎每个扫描点的噪点都会被当成极窄峰识别出来。还有可能是ppm设太大把同位素峰、加合峰都当成独立候选了。提示判断峰数量是否合理没有绝对标准。一般做正离子模式的血浆血清样本一个样本检测到5000~20000个峰是常见的。如果只有几百个大概率是参数太严了如果动不动几十万个那多半是模型对噪声过拟合了。5.2 内存溢出或跑批速度慢xcms最让人头疼的另一个点就是当你有一百多个大体积.mzML文件时跑起来非常吃内存。我的建议是数据量大的时候不要一次性把所有文件都读进R。先对每个文件单独做峰检测再把检测结果合并。xcms的findChromPeaks本身支持传入一个文件列表并逐个处理但要注意合并时的分块策略。还有一招是把数据切成多个时间窗口分段处理最后再用joinWindows合并这招对保留时间很长的梯度特别有效。另外能开多线程就开多线程。我在Linux服务器上跑register(BPPARAM MulticoreParam(4))直接提速三倍以上。5.3 保留时间校正后反而更差了有一次我拿一批跨了三个月、分三批采集的样本跑数据发现adjustRtime做完之后部分峰的对齐效果反而比校正前更差。后来排查发现是因为样本间空白批次太多空白样本里的峰很稀疏整个校正过程被这些空白样本严重带偏。这时候有两个思路一是把空白样本从校正流程中暂时剔除只用真实样本做对齐之后再单独补上空白样本的峰二是用ObiwarpParam时把distFun从默认的cor换成cor_opt或cov后两者在稀疏数据下更稳健。这个方法调了之后峰表的完整度明显好了很多。5.4 加合物和同位素峰导致特征冗余xcms输出的峰表里面同一个代谢物往往对应多个特征。比如一个化合物有MH峰还有MNa峰还有MH的13C同位素峰。如果不做去冗余后面做差异代谢物筛选时会看到一堆高度相关的重复特征统计上容易造成假阳性。这边我的经验是按相关性聚类的思路来处理先根据峰表计算特征间的相关矩阵然后用一个联通分支算法把相关性大于0.8的特征聚成一组每组里保留信号最强或最稳定的那个作为代表特征。更精细的办法是用CAMERA另一个包来识别加合物和同位素峰它能把同位素峰、加合峰归到一个“feature group”里。不过注意CAMERA本身的输出格式和xcms新版本略有出入用之前最好查一下兼容性。5.5 新版本xcms和老代码的兼容性xcms的API在过去十年里大改过好几次。你如果是在网上搜的教程看到有人用xcmsSet、group、retcor这些老函数千万别直接复制。新版本3.x推荐的写法是用findChromPeaks、groupChromPeaks、adjustRtime这套函数。如果你拿到的是老代码可以粗略理解成new(xcmsSet)→readMSDatafindPeaks.centWave()→findChromPeaks()group()→groupChromPeaks()retcor()→adjustRtime()getPeaklist()→featureDefinitions()featureValues()为了省事我现在直接统一用findChromPeaks、groupChromPeaks、adjustRtime这套新接口跑完再用exportMetaboAnalyst之类的方式导出可以和下游工具无缝衔接。6. 一个可参考的完整代码骨架最后我把一套在实际项目里跑通的完整代码骨架贴出来你可以直接照着改成自己的样式。注意这套代码只是参考参数务必按你自己的数据特性来调。library(xcms) library(MSnbase) # 1. 读入数据 raw_files - list.files(mzML_dir, pattern .mzML, full.names TRUE) raw_data - readMSData(raw_files, mode onDisk) # 2. 设置分组信息假设文件名condition_rep.mzML sample_info - data.frame( sample basename(raw_files), group sub(_.*, , basename(raw_files)) ) sampleData(raw_data)$group - sample_info$group # 3. 峰检测 cwp - CentWaveParam( peakwidth c(5, 20), snthresh 10, ppm 15, noise 1000, fitgauss TRUE, verboseColumns TRUE ) xdata - findChromPeaks(raw_data, param cwp) # 4. 初始分组用于校正 pdp - PeakDensityParam( sampleGroups sample_info$group, bw 5, minFraction 0.5, binSize 0.025 ) xdata - groupChromPeaks(xdata, param pdp) # 5. 保留时间校正 xdata - adjustRtime(xdata, param ObiwarpParam()) # 6. 最终分组 xdata - groupChromPeaks(xdata, param pdp) # 7. 填充缺失峰 xdata - fillChromPeaks(xdata, param FillChromPeaksParam()) # 8. 提取峰表和特征矩阵 peak_table - chromPeaks(xdata) feature_matrix - featureValues( xdata, method maxint, value into, filled TRUE ) # 9. 保存 write.csv(feature_matrix, xcms_feature_matrix_final.csv, quote FALSE)7. 我的实操体会与扩展建议用了这么多年现在回头再看xcms我依然觉得它“正中代谢组学的眉心”这个说法并不过分。它没有做什么天花乱坠的事情但把数据处理中最基础、最枯燥、最繁琐的工作用一套开放且可继承的框架固化了下来。尤其是当你需要批量处理上百个样本或者需要完全重现一套数据处理流程的时候那种“一切都清清楚楚每一步都有日志”的踏实感是任何商业软件都给不了的。如果你现在刚开始学代谢组学数据处理我建议不要跳过xcms直接去点各种在线平台。花两周时间把xcms的基本工作流跑通你就能深刻理解峰检测、保留时间校正、峰对齐这些名词的实际意义。等回过头再去看其他软件你会发现它们说的很多“改进”和“优势”底层逻辑其实都是在和xcms做对比。最后提醒一句跑数据前一定要先看好仪器状态同一天连续采集的样本质谱漂移小参数就好调跨很多天采集的样本参数就得留出足够裕量。数据处理永远是为数据采集服务的设备端做不好算法再优化也白搭。