代谢组学数据清洗这件事说难不难说简单也真不简单。我最早接触这批需求是在处理一批非靶向代谢组学的注释结果时手里攥着几百上千条化合物记录每条都带着 Formula分子式和 INCHIKEY但下游分析要用的 HMDB ID 和 Common Name 却是空的。手工一条条去查不现实。写脚本硬匹配又会撞上一堆坑——分子式对上了但结构不对、INCHIKEY 前 14 位相同却分属不同条目、同一个化合物在不同数据库里叫法五花八门。这篇就围绕用 Formula 值匹配 INCHIKEY再拿到对应 HMDB ID 和 Common Name这条链路把我踩过的坑、验证过的思路和能直接抄的代码完整拆一遍适合做代谢组学、天然产物分析、化合物数据库整合的朋友参考新手能照着跑通老手也能在细节里找到几个省时间的点。1. 先搞清楚这条匹配链路到底在解决什么问题1.1 Formula、INCHIKEY、HMDB ID 三者各自的角色很多人一上来就想写 SQL 或者 pandas 的 merge结果发现字段根本对不上。问题出在没理清这三个字段的语义层级。Formula分子式是最粗的一层它只描述原子组成比如 C6H12O6对应的可能是葡萄糖、果糖、半乳糖等一大堆同分异构体。INCHIKEY是 InChI 的哈希压缩形式标准长度 27 位前 14 位编码的是分子骨架连接性也就是这个分子长什么样后 10 位编码的是立体化学、同位素、电荷等细节最后一位是校验位。HMDB ID则是人类代谢组数据库里的唯一主键形如 HMDB0000122它绑定的是一个具体的数据库条目附带 Common Name、分类、生物样本浓度等一堆元数据。所以这条链路的本质是用分子式做粗筛用 INCHIKEY 做精确定位最终落到 HMDB 条目上取名字。分子式只是入口INCHIKEY 才是真正的身份凭证。我见过不少人直接拿 Formula 去 HMDB 里模糊搜结果一个 C6H12O6 返回几十条根本没法自动判断该选哪个——这就是跳过 INCHIKEY 的代价。1.2 为什么不能只用 Formula 直接查 Common Name举个实际场景。你有一批数据Formula 列写着 C10H16N5O13P3这是 ATP 的分子式。但如果你只拿这个分子式去匹配数据库里可能同时存在 ATP、dATP 的近似条目、各种磷酸化衍生物甚至一些分子式相同但结构完全不同的代谢物。Common Name 一旦匹配错下游的通路富集分析、生物标志物筛选全部跟着错而且这种错误极其隐蔽——名字看起来都像那么回事。INCHIKEY 的价值就在这里它是结构决定的同一个 INCHIKEY 前 14 位基本锁定了分子的连接方式能有效排除同分异构体的干扰。所以正确的顺序永远是Formula 辅助校验 INCHIKEY 主键匹配而不是反过来。我在实际项目里会把 Formula 当作一道防呆检查如果某条记录的 INCHIKEY 匹配到了 HMDB 条目但该条目的分子式和输入的 Formula 对不上那这条匹配就要打上可疑标记人工复核。这一步能拦下不少因为 INCHIKEY 录入错误或数据库版本差异导致的误配。1.3 这条链路典型的应用场景说几个我实际遇到过的用法。第一类是非靶向代谢组学注释结果的下游整理质谱软件输出的化合物表里往往只有 Formula 和 INCHIKEY缺 HMDB ID而后续做通路分析、网络分析又必须要 HMDB ID 做索引。第二类是多数据库交叉验证你从 PubChem 或 ChemSpider 拿到一批化合物想映射到 HMDB 体系里看它们的生物学角色。第三类是自建化合物库的标准化实验室积累的化合物清单字段混乱需要统一补全标准 ID 和通用名。这三类场景的共同点是数据量大、人工不可行、对准确性要求高。所以下面讲的方案不是能跑就行而是要考虑可复现、可审计、出错能定位。2. 数据源怎么选本地库、在线接口还是现成映射表2.1 三种数据获取路径的取舍做 INCHIKEY 到 HMDB 的映射数据从哪来直接决定了方案的稳定性和速度。我梳理了三条常见路径各有各的脾气路径数据来源优点缺点适用场景本地 HMDB 结构文件HMDB 官方提供的 XML/SDF 全量下载离线、可批量、字段完整文件大、更新需重新下载、解析成本高大批量、重复性任务在线化合物接口各类化合物数据库的查询 API实时、无需维护本地库有速率限制、网络依赖、字段不全小批量、临时查询现成映射表社区整理的 INCHIKEY-HMDB 对照表开箱即用、轻量覆盖不全、版本滞后、来源需甄别快速验证、补充我的建议是以本地 HMDB 结构文件为主在线接口为辅。原因很实在代谢组学项目动辄几千条记录走在线接口光是速率限制就能把你卡到怀疑人生而且一旦网络抖动整个流程就得重跑。本地库虽然初次准备麻烦点但一次建好后面所有项目都能复用长期看时间成本低得多。2.2 本地 HMDB 结构文件的字段结构HMDB 官方提供的结构文件里每条代谢物记录通常包含这些关键字段HMDB ID、Common Name、Chemical Formula、Average Molecular Weight、SMILES、InChI、InChIKey、Synonyms 等。解析的时候重点抓四个accession即 HMDB ID、nameCommon Name、chemical_formula、inchikey。这里有个容易忽略的细节同一个 HMDB 条目可能对应多个 INCHIKEY尤其是带立体异构的化合物数据库里会记录不同立体形式的 INCHIKEY。反过来同一个 INCHIKEY 前 14 位可能对应多个 HMDB 条目不同立体化学或不同盐形式。所以建索引的时候我通常建两层一层用完整 27 位 INCHIKEY 精确匹配一层用前 14 位骨架匹配做兜底。两层配合既能保证精度又能在完整键对不上时给出候选。2.3 在线接口作为补充时的注意事项在线接口不是不能用但要清楚它的边界。第一速率限制必须做请求间隔和重试机制我一般设 0.2 到 0.5 秒的间隔遇到限流就指数退避。第二返回字段可能不全有些接口只给 HMDB ID 不给 Common Name那就得二次查询。第三结果要缓存同一个 INCHIKEY 查过一次就落本地避免重复请求。这三点做到了在线接口作为补充手段是够用的。提示无论用哪种数据源都要记录数据版本和获取时间。代谢组学数据库会更新今天匹配上的条目明天可能被合并或修订没有版本记录的话几个月后你根本说不清当时的结果是怎么来的。3. 匹配逻辑的核心设计从精确到模糊的分层策略3.1 第一层完整 INCHIKEY 精确匹配这是最理想的情况。输入的 INCHIKEY 是完整 27 位本地库里也有完全一致的记录直接一对一映射拿到 HMDB ID 和 Common Name。这一步的代码逻辑非常简单本质就是一个字典查找或者 DataFrame 的 merge。但实际数据里完整 27 位能精确命中的比例往往没有想象中高。原因有几个输入数据的 INCHIKEY 可能被截断、可能带了多余的空格或大小写不一致、可能是旧版本算法生成的。所以匹配前必须做标准化清洗去首尾空格、统一转大写、检查长度是否为 27、校验最后一位校验位。我见过太多因为一个尾随空格导致整批匹配失败的案例这种坑完全可以在预处理阶段消灭。3.2 第二层前 14 位骨架匹配做兜底当完整键对不上时退而求其次用前 14 位匹配。前 14 位锁定的是分子连接性能保证骨架相同但立体化学可能不同。这时候匹配结果可能是一对多需要额外逻辑来裁决。我的处理方式是前 14 位匹配到的所有候选再用 Formula 做一次交叉验证。如果候选条目的分子式和输入 Formula 一致保留不一致的剔除。如果交叉验证后还剩多个候选就标记为多候选待定输出所有候选的 HMDB ID 和 Common Name 供人工判断而不是武断地取第一个。这一点很重要——自动流程最怕的就是看起来成功了但其实是错的宁可标记待定也不要制造假阳性。3.3 第三层Formula 反向校验与异常标记前面说过 Formula 是防呆检查具体怎么落地我的做法是在匹配完成后对每一条成功映射的记录做一次反向核对拿匹配到的 HMDB 条目的分子式和输入的 Formula 做字符串比对。这里要注意分子式的书写规范问题——元素顺序、大小写、括号使用都可能不一致比如 C6H12O6 和 H12C6O6 在字符串层面不等但化学上等价。所以比对前要做分子式规范化解析成元素计数字典再按标准顺序通常是 C、H 然后其他元素字母序重新拼接。规范化之后比对能消除大部分书写差异带来的假阴性。对于比对不上的记录打上formula_mismatch标记单独输出一个异常清单。这个清单在实际项目里非常有用往往是数据录入错误或者数据库版本差异的线索。3.4 匹配优先级的完整决策流程把三层逻辑串起来完整的决策流程是这样的输入记录先做 INCHIKEY 标准化去空格、转大写、长度校验。用完整 27 位键查本地索引命中则直接取 HMDB ID 和 Common Name标记match_level exact。未命中则用前 14 位查骨架索引得到候选集。候选集用规范化后的 Formula 过滤剩一条则标记match_level skeleton剩多条则标记match_level ambiguous并输出全部候选。所有命中记录做 Formula 反向校验不一致的打formula_mismatch标记。完全未命中的记录标记match_level none进入人工或在线接口补充流程。这个流程的好处是每一步的匹配质量都有明确标记后续做统计分析时可以按 match_level 分层处理比如只保留 exact 和 skeleton 级别的结果做严格分析ambiguous 的单独讨论。4. 用 Python 把整条链路跑通可复现的实操代码4.1 环境准备与依赖这套流程用 Python 就够了核心依赖不多pandas处理表格lxml或xml.etree解析 HMDB 的 XML 文件requests用于可选的在线补充。装依赖这一步没什么好省的直接上pip install pandas lxml requests如果你的 HMDB 结构文件是 SDF 格式还需要rdkit来解析但纯做 INCHIKEY 映射的话XML 版本字段更规整我更推荐用 XML。SDF 解析虽然灵活但字段提取要写不少样板代码对只想拿 ID 和名字的场景属于杀鸡用牛刀。4.2 解析 HMDB 结构文件建本地索引假设你拿到的是 HMDB 的代谢物 XML 文件结构大致是每条metabolite节点下挂着accession、name、chemical_formula、inchikey等子节点。解析的核心是把这些字段抽出来建成两个字典索引import xml.etree.ElementTree as ET import pandas as pd def build_hmdb_index(xml_path): tree ET.parse(xml_path) root tree.getroot() records [] for metab in root.iter(metabolite): accession metab.findtext(accession) name metab.findtext(name) formula metab.findtext(chemical_formula) inchikey metab.findtext(inchikey) if inchikey: inchikey inchikey.strip().upper() records.append({ hmdb_id: accession, common_name: name, formula: formula, inchikey: inchikey }) df pd.DataFrame(records) # 完整键索引 exact_index {} for _, row in df.iterrows(): if row[inchikey] and len(row[inchikey]) 27: exact_index.setdefault(row[inchikey], []).append(row) # 骨架索引前14位 skeleton_index {} for _, row in df.iterrows(): if row[inchikey] and len(row[inchikey]) 14: skeleton_index.setdefault(row[inchikey][:14], []).append(row) return df, exact_index, skeleton_index这里用setdefault配合列表是因为一个键可能对应多条记录必须保留全部候选而不是覆盖。这一点在骨架索引上尤其关键前 14 位一对多是常态。4.3 分子式规范化函数前面强调过分子式比对前要规范化具体实现是把分子式解析成元素计数字典再重排。这里要处理括号和嵌套的情况比如 Ca(OH)2 这种。写一个能处理常见情况的解析器import re from collections import defaultdict def normalize_formula(formula): if not formula: return None formula formula.strip() # 处理括号展开 def expand(f): # 简单括号处理支持一层嵌套 pattern re.compile(r\(([^()])\)(\d*)) while pattern.search(f): f pattern.sub(lambda m: m.group(1) * int(m.group(2) or 1), f) return f expanded expand(formula) tokens re.findall(r([A-Z][a-z]?)(\d*), expanded) counts defaultdict(int) for elem, num in tokens: counts[elem] int(num) if num else 1 # 标准顺序C、H 优先其余字母序 ordered [] for e in [C, H]: if e in counts: ordered.append(f{e}{counts[e] if counts[e] 1 else }) for e in sorted(k for k in counts if k not in (C, H)): ordered.append(f{e}{counts[e] if counts[e] 1 else }) return .join(ordered)这个函数不是万能的遇到非常复杂的嵌套括号或水合物写法可能出错但对代谢组学里绝大多数分子式够用了。实际使用时如果规范化失败返回 None 或抛异常就把该条记录标记为formula_parse_error不要让它静默通过。4.4 主匹配函数与结果输出把索引和规范化函数组合起来写主匹配逻辑def match_records(input_df, exact_index, skeleton_index): results [] for _, row in input_df.iterrows(): raw_key str(row.get(INCHIKEY, )).strip().upper() input_formula normalize_formula(row.get(Formula, )) result { input_inchikey: raw_key, input_formula: row.get(Formula, ), hmdb_id: None, common_name: None, match_level: none, formula_mismatch: False } if len(raw_key) 27 and raw_key in exact_index: cand exact_index[raw_key][0] result.update({ hmdb_id: cand[hmdb_id], common_name: cand[common_name], match_level: exact }) if input_formula and normalize_formula(cand[formula]) ! input_formula: result[formula_mismatch] True elif len(raw_key) 14 and raw_key[:14] in skeleton_index: cands skeleton_index[raw_key[:14]] if input_formula: filtered [c for c in cands if normalize_formula(c[formula]) input_formula] else: filtered cands if len(filtered) 1: result.update({ hmdb_id: filtered[0][hmdb_id], common_name: filtered[0][common_name], match_level: skeleton }) elif len(filtered) 1: result[match_level] ambiguous result[candidates] ;.join( f{c[hmdb_id]}|{c[common_name]} for c in filtered) results.append(result) return pd.DataFrame(results)输出结果里match_level和formula_mismatch两个字段是后续人工复核的关键抓手。我一般会把结果按 match_level 分组导出成多个 sheetexact 和 skeleton 的直接用ambiguous 和 none 的单独处理。4.5 结果验证抽样回查与统计代码跑完不代表结束必须做验证。我的习惯是按 match_level 分层抽样每层随机抽 20 到 30 条拿 HMDB ID 去数据库里回查确认 Common Name 和 Formula 都对得上。同时统计各层的占比如果 exact 占比异常低比如低于 30%就要怀疑输入数据的 INCHIKEY 质量或者本地库版本是否匹配。还有一个统计值得看formula_mismatch的比例。如果这个比例超过 5%说明要么输入 Formula 有系统性问题比如格式不统一要么本地库和输入数据来自不同版本需要排查。这些统计数字不是走形式它们是判断整批数据可信度的直接依据。5. 那些文档里不会写的坑我踩过的五个典型问题5.1 INCHIKEY 大小写与隐藏字符这是最高频的坑没有之一。从 Excel 复制出来的 INCHIKEY 经常带不可见的尾随空格或者制表符肉眼完全看不出来但字符串比对就是失败。更隐蔽的是大小写混用有些软件输出小写有些输出大写而字典查找是大小写敏感的。解决办法只有一个匹配前统一strip().upper()并且用正则检查是否只含合法字符。我现在的习惯是在预处理阶段就打印出所有长度不等于 27 的记录逐条看往往能揪出一批脏数据。5.2 同分异构体导致的骨架匹配歧义前 14 位匹配一对多是设计内的行为但歧义怎么裁决是个真问题。我遇到过一批糖类化合物骨架相同、Formula 相同但立体构型不同前 14 位完全一样Formula 也过滤不掉。这种情况下任何自动裁决都是猜正确做法是标记 ambiguous 并输出全部候选让懂化学的人来判断。强行取第一个的后果是引入系统性错误而且因为名字看起来合理极难被发现。5.3 分子式书写顺序不一致引发的假阴性C6H12O6 和 H12C6O6 化学上等价字符串比对却不等。如果不做规范化反向校验会误报大量 formula_mismatch把本来正确的匹配标记成可疑。这个坑的教训是任何涉及化学式比对的环节都必须先规范化。规范化函数要能处理元素顺序、括号、数字省略单个原子不写 1这几种常见差异。5.4 数据库版本差异导致的 ID 失效HMDB 会定期更新条目可能被合并、拆分或废弃。你三个月前建的本地索引今天拿来匹配新数据可能就有条目对不上了。更麻烦的是有些旧 ID 在新版本里已经不存在但输入数据里还带着。我的应对是在索引里保留条目的状态字段如果源文件有的话匹配时如果命中已废弃条目标记出来提醒更新。同时记录本地库的版本号和下载日期写进结果文件的元数据里。5.5 一对多映射时的输出设计一个输入 INCHIKEY 可能对应多个 HMDB 条目不同盐形式、不同立体化学这时候输出格式的设计很关键。如果只输出一个字段多候选信息就丢了如果拆成多行下游统计又会重复计数。我的做法是主表保持一行一条输入记录多候选信息压缩到一个字段里用分隔符拼接需要展开时再单独解析。这样既保证主表行数稳定又不丢信息。6. 匹配质量的分层处理与人工复核策略6.1 按 match_level 分层导出与优先级跑完匹配后我从不把结果一股脑丢给下游而是按 match_level 分层导出。exact 级别的结果可信度最高可以直接进入分析流程skeleton 级别的次之建议抽样验证后使用ambiguous 和 none 的必须人工介入。这种分层不是形式主义它让整个流程的可审计性大大提升——出了问题能快速定位是哪一层出的错。具体操作上我会生成一个汇总统计表列出各 match_level 的记录数、占比以及 formula_mismatch 的数量。这个表放在结果文件的第一页任何人拿到数据先看这个表就能对整体质量有个判断。6.2 ambiguous 候选的人工裁决辅助对于 ambiguous 的记录光给候选列表还不够得给裁决提供依据。我会在候选信息里额外附上每个候选的分子量、SMILES如果有的话方便复核人员快速判断。如果候选数量太多比如超过 5 个说明前 14 位骨架太通用这时候建议改用完整键或者结合其他字段如保留时间、质荷比来缩小范围。人工裁决的结果要回写到主表并且记录裁决人和裁决依据。这一步在团队协作里尤其重要否则过几个月没人记得当时为什么选了某个候选。6.3 未匹配记录的补充查询流程match_level 为 none 的记录先别急着放弃。我的流程是先用在线接口按 INCHIKEY 查一遍如果在线也查不到再用 Formula 加分子量做模糊查询看是否有近似条目。在线查询的结果同样要缓存和记录来源。如果所有途径都失败就标记为unresolved在报告里单独说明数量而不是悄悄丢掉——未匹配记录的比例本身就是一个数据质量指标。7. 把这套流程用顺之后的几点体会这套 Formula 加 INCHIKEY 匹配 HMDB 的流程我从最早的手忙脚乱到现在能稳定批量跑中间最大的转变不是代码写得多漂亮而是对数据质量的敬畏。刚开始我总想着一步到位全自动结果被各种脏数据教做人。后来想明白了自动化的目标不是消灭人工而是把人工集中在真正需要判断的地方。exact 匹配交给机器ambiguous 和异常交给人工各司其职整体效率和准确性反而都上去了。另外一个体会是索引和规范化这两个基础环节值得多花时间打磨。索引建得好匹配逻辑就简单规范化做得全假阴假阳就少。很多人急着写匹配主逻辑结果基础不牢后面全在打补丁。我现在的习惯是先把数据源摸清楚、把字段清洗规则定死、把索引结构设计好再动手写匹配前期多花的两小时能省后面两天的调试。最后分享一个小技巧每次跑完批量匹配我都会把结果和上一批的统计数字对比一下。如果 exact 占比突然下降或者 formula_mismatch 突然上升多半是输入数据或者本地库出了问题这时候停下来查比闷头往下跑划算得多。数据清洗这活儿慢就是快。