4、SMILES字符串解析:如何读取和解析SMILES
拿到SMILES之后,接下来就是怎么读、怎么拆的问题了。我个人习惯把SMILES解析叫做“拆积木”,因为本质上它就是把一个字符串里的化学信息还原成原子和键的过程。
4.1 SMILES最基本的读取思路
其实SMILES解析没有你想象中那么玄乎。我当年第一次写解析器时,觉得这是要懂编译原理的活儿,后来发现核心逻辑就是个状态机配上栈。
说白了,你需要做到这几步:
- 从左往右一个个字符扫过去;
- 遇到字母就是原子;
- 遇到数字就是环编号;
- 遇到括号就是分支;
- 遇到方括号就是带方括号的原子(比如[Fe]、[OH-]);
- 遇到等号、井号、冒号就是键类型。
嗯,这里要注意,SMILES里默认的单键是不写出来的。C-O这个写法虽然没写“-”,但你知道中间是一定有单键的。刚入门的同学经常在这块脑补过度,其实SMILES默认规则就是你“看见相邻就算单键”,除非你明确写了别的符号。
4.2 解析过程中的几个核心技巧
我的做法是先做一个词法分析,也就是先把SMILES文本拆成一个一个的tokens。这个办法我在早期项目里吃了不少苦头才学会——直接用字符串操作去处理分支和环,很快就炸了。
比如下面这段代码,就是把SMILES拆成基础token的骨架逻辑:
def tokenize(smiles: str):
tokens = []
i = 0
while i < len(smiles):
ch = smiles[i]
if ch.isalpha(): # 原子符号
if i+1 < len(smiles) and smiles[i+1].islower():
tokens.append(smiles[i:i+2])
i += 2
else:
tokens.append(ch)
i += 1
elif ch in "()[]": # 括号、分支
tokens.append(ch)
i += 1
elif ch.isdigit(): # 环编号
tokens.append(ch)
i += 1
elif ch in "=#:/-\\@": # 键或手性符号
tokens.append(ch)
i += 1
else:
i += 1
return tokens
这段代码我反复改过。有些极端案例比如[C@@H]这种手性标记,字符长度会超出两个,但基础的tokenizer足够覆盖常见的SMILES。
解析的第二步是建分子图。你拿到tokens以后,需要搭配一个栈(Python的list就行)来跟踪当前的原子节点。遇到左括号(就是分支入口)我就push当前节点,遇到右括号就pop回去。
4.3 我踩过的几个解析坑(避坑指南)
第一个坑:原子符号的大小写。
我曾经有一回,用smiles[i].isupper()来判断原子符号的结尾,结果Co(钴)这种双字母元素全给拆成了C和o。Co变成了一个碳加上一个氧——这化学味道都不对了。
解决方案:往后多看一眼,if smiles[i].isupper() and smiles[i+1].islower()才走双字母分支。
第二个坑:环的编号对应。
SMILES里环编号上首次出现是开环,第二次出现是闭环。但编号可以打乱顺序写,比如C1CC2C1C2是可以的。你在解析时要维护一个ring_map,第一次遇到数字的时候把当前节点存起来,第二次遇到同一个数字时就把两个节点连一条键。
我这里用一个简化的思路:
# 伪代码
ring_map = {}
for token in tokens:
if token.isdigit():
if token not in ring_map:
ring_map[token] = current_atom
else:
connect_ring(current_atom, ring_map[token])
del ring_map[token]
第三个坑:分支里套分支。
你看这个SMILES:C(C(C)C)C,三层的嵌套。如果你解析时没用好栈,就会跑到某个括号里就回不来了。我写过一版没有用栈的代码,结果挂了十几个公司内部测试用例,全是这种嵌套分支。
4.4 完整的SMILES解析流程
整体流程可以用下面这张图概括。你先做预处理,然后词法扫描,再用栈+原子链表建分子图,最后检查一下环有没有闭合、价键有没有超限。
R
4.5 主流Python工具是怎么解析的
实际项目里,我一般不自己写解析器,而是依赖于几个成熟的库。因为自己做容易漏掉手性标记和同位素这些细节。
下面是三个最常见的选项,我用表格列一下差异:
| 库 | 解析方式 | 个人评价 |
|---|---|---|
| RDKit | C++底层调boost::graph,Python封装 | 稳、快,但安装时依赖不少 |
| Open Babel | 纯C++ + SWIG调用 | 历史久,生态杂,实测SMILES兼容率99%+ |
| 自写解析器 | 纯Python手工tokenizer | 教学和轻量场景可以用,但不建议上生产 |
RDKit是我现在的主力。它的解析器经历了十几年开源社区的考验。我之前在一个药物筛选项目里跑了1300万条SMILES,RDKit只报错了不到10条,而且全是SMILES本身非法。
4.6 解析后的检验与容错
解析完SMILES并不代表万事大吉。我通常会再加四道关:
- 检查价电子是否匹配——碳连了五个氢肯定有问题;
- 检查环有没有闭合——如果有只开不闭的数字,多半SMILES是半残的;
- 检查括号是否匹配——左括号数量必须等于右括号;
- 检查手性中心——像@和@@这些符号不能重复也不能漏写。
Chem.MolFromSmiles() 返回了None,而我之前写的检查流程直接把它当空值处理掉,结果导致了分子信息遗漏。后来我加了一条规则:解析失败时记录日志并保留原文,不做静默过滤。
4.7 解析结果的常见误区
解析完之后,SMILES字符串里的原子顺序跟你预期的不一定一样。RDKit的SMILES输出的原子顺序是按原子在遍历树时的访问次序来的,不是按字符串长度或分子量。这个我刚开始做分子指纹的时候迷惑了好一阵子,因为排好序的分子图和原始SMILES看起来完全不是一回事,但化学结构完全正确。
所以,你读入SMILES后的第一个操作,最好是用Chem.CanonSmiles(smiles)做一个规范化。这样可以保证同样的分子只有一种表示方式。具体用法是:
from rdkit import Chem
canon_smi = Chem.MolToSmiles(Chem.MolFromSmiles("C1=CC=CC=C1"))
print(canon_smi) # 输出: c1ccccc1
你看,同一个苯环,带数字的C1=C风格进去,出来的就是小写的c1ccccc1规范化形式。这种统一格式在后续做序列长度统计时至关重要——你不想因为SMILES的写法差异就多算或少算一个字符吧。
做可视化教程时,我最常用的一句代码就是先规范、再统计、最后画图。你想想看,要是没有这层解析和规范,后面所有关于SMILES长度的结论都是错的。