分子中原子编号的匹配|Jerkwin

2018-01-25 06:43:45

以前的时候我曾简单介绍过有关分子匹配的问题分子匹配叠合及两种冰的结构异同, 也曾示例过如何借助GaussView来实现分子片段分析. 但所用的方法总觉得不很方便, 所以就想着自己能不能实现一个类似功能的简单脚本或在线工具. 但正如以前所说, 就一般情况而言, 这并不是一件简单的问题. 不过针对我们需要的功能, 还是有一些算法可以实现的, 主要是图匹配算法.

图匹配算法

如果你也需要了解这方面的知识, 请先阅读下面的相关资料.

有没有什么算法可以确定两图是否同构？
如何通过邻接矩阵判断拓扑结构是等价的？
Reading List: Graph Isomorphism
One of These Things is Not Like The Others
Revised algorithm (with a lot cleaner source code) that covers regular and weird graphs
Ullman’s Subgraph Isomorphism Algorithm
VF2 algorithm steps with example
图的同构和 Girth 问题
（子）图同构算法VF2实现（2）——源代码与具体实现
Xiao-xiang Zou; 图同构中的一类顶点细分方法; Journal of Software 18(2)：213, 2007; 10.1360/jos180213
H. L. Morgan; The Generation Of A Unique Machine Description For Chemical Structures-a Technique Developed At Chemical Abstracts Service.; J. Chem. Doc. 5(2)：107-113, 1965; 10.1021/c160017a018
Brendan D. Mckay, Adolfo Piperno; Practical Graph Isomorphism, Ii; Journal of Symbolic Computation 60：94-112, 2014; 10.1016/j.jsc.2013.09.003
宋婷婷; 工程化学数据库中分子结构信息的处理
唐玉焕, 林克江, 尤启冬; 基于2D分子指纹的分子相似性方法在虚拟筛选中的应用
余靖, 韩玉; Ursi：高效的子图同构查询算法
王利莎, 袁身刚, 欧阳政, 郑崇直; 目标化合物析分系统中的重要算法

总结一下, 目前常用的图匹配算法包括, Ullman, VF2, NAUTY, 效率依次递增, 但算法复杂度也依次递增.

可以参考的实现代码如下:

js的Ullman: subgraph-isomorphism
python的VF2: isomorphvf2.py
c的NAUTY: nauty and Traces或旧版本建议使用2.0版本. 最新版本功能太多, 太复杂, 不适合集成到自己的代码中

在匹配时, 还可以先使用化学信息学中常用的Morgan算法做粗略筛选, 但它不能保证匹配结果一定正确.

程序实现

当两个相同的分子具有不同的构型和原子编号顺序时, 如果我们要对它们的原子编号进行匹配, 在实现的时候不仅要考虑图匹配的作法, 还要考虑每个原子的元素种类和化学环境, 因此, 相比上面所讨论的图匹配, 我们还要适当考虑问题的特殊性. 通过提前对原子使用元素与成键原子类型进行匹配筛选, 可以大大减少后面图匹配的可能性, 提高匹配速度.

为了方便做成在线工具, 我选择了使用已有js代码的Ullman算法. 但事实证明, 这个算法耗时随原子数指数增大. 当分子中的原子数超过20之后, 算法运行时间过长, 基本没有意义, 因此这个算法的实用性不好. 为此, 我在程序中增加了手动添加”锚点”原子的功能, 也就是手动指定两个分子中必须匹配的原子编号. 利用锚点原子, 可以大大加快匹配速度. 很显然, 锚点原子越多, 匹配越迅速. 当锚点原子数与待匹配原子数相同时, 就相当于手动指定每个匹配原子了. 根据我的测试, 在不使用锚点原子的情况下, 30个原子的匹配已经是算法极限了, 但使用锚点原子后, 处理50个原子的匹配不成问题. 至于更多原子的情况, 还有待测试.

使用简介

我将这个工具命名为matchmol. 它是一个在线工具, 在浏览器中运行. 目前仅支持读入GaussView格式的PDB文件, 其中必须包含元素符号列和成键列表, 因为需要使用这些信息进行匹配.

加载两个PDB文件, 其分子相同, 但具有不同的原子编号顺序. 左方为参考分子, 右方为待匹配分子, 中间的文本框中列出每个分子的原子编号, 参考分子的编号顺序不可更改, 待匹配分子的编号可以更改.

可使用三种方法将待匹配分子的原子编号顺序调整至与参考分子一致:

(全)自动模式: 点击Auto Match进行自动匹配. 如果找到匹配结果, 则会在右方文本框中列出待匹配分子中与参考分子对应的编号顺序. 点击Application, 则会根据列出的编号顺序调整待匹配分子的原子编号, 方便对比. 如果点击Save PDB File则会将调整顺序后的待匹配分子的PDB文件保存为conf.pdb. 注意: 全自动模式能处理的分子其原子数在20以下, 否则运行时间过长, 容易失败.
锚点模式(半自动模式): 如果待匹配原子数过多, 自动模式无法在短时间内完成匹配, 可手动指定锚点原子, 加速匹配. 指定方法是在待匹配分子原子编号文本框中使用带有*的编号. 举例来说, 如果参考分子中的3号原子对应于待匹配分子的7号原子, 那么将待匹配分子原子编号文本框中的3改为7*, 即表示3号原子锚定为7号原子. 接下来的操作与自动模式相同.
手动模式: 如果你手动指定每个原子的锚点原子, 那就无须在编号后面加*了. 注意, 这种情况下必须指定所有原子的锚点原子.

应用示例

简单分子

以下为同一分子不同原子编号顺序下的PDB文件, 用于说明支持的PDB文件的格式.

mol1.pdb
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33	TITLE Molecule Name REMARK 1 File created by GaussView 5.0.9 HETATM 1 C 0 -2.958 -0.067 -2.907 C HETATM 2 H 0 -2.409 -0.788 -2.338 H HETATM 3 C 0 -3.408 -0.697 -4.233 C HETATM 4 O 0 -3.739 -1.861 -4.248 O HETATM 5 N 0 -3.427 0.073 -5.380 N HETATM 6 C 0 -3.804 -0.532 -6.661 C HETATM 7 H 0 -4.713 -1.084 -6.541 H HETATM 8 H 0 -3.948 0.238 -7.390 H HETATM 9 C 0 -3.090 1.494 -5.474 C HETATM 10 H 0 -3.824 1.996 -6.069 H HETATM 11 H 0 -3.073 1.922 -4.493 H HETATM 12 H 0 -3.027 -1.191 -6.987 H HETATM 13 H 0 -2.127 1.603 -5.928 H HETATM 14 H 0 -2.334 0.779 -3.108 H HETATM 15 H 0 -3.817 0.247 -2.352 H END CONECT 1 3 2 14 15 CONECT 2 1 CONECT 3 1 4 5 CONECT 4 3 CONECT 5 3 6 9 CONECT 6 5 7 8 12 CONECT 7 6 CONECT 8 6 CONECT 9 5 10 11 13 CONECT 10 9 CONECT 11 9 CONECT 12 6 CONECT 13 9 CONECT 14 1 CONECT 15 1

mol1.pdb

TITLE      Molecule Name
REMARK   1 File created by GaussView 5.0.9
HETATM    1  C           0      -2.958  -0.067  -2.907                       C
HETATM    2  H           0      -2.409  -0.788  -2.338                       H
HETATM    3  C           0      -3.408  -0.697  -4.233                       C
HETATM    4  O           0      -3.739  -1.861  -4.248                       O
HETATM    5  N           0      -3.427   0.073  -5.380                       N
HETATM    6  C           0      -3.804  -0.532  -6.661                       C
HETATM    7  H           0      -4.713  -1.084  -6.541                       H
HETATM    8  H           0      -3.948   0.238  -7.390                       H
HETATM    9  C           0      -3.090   1.494  -5.474                       C
HETATM   10  H           0      -3.824   1.996  -6.069                       H
HETATM   11  H           0      -3.073   1.922  -4.493                       H
HETATM   12  H           0      -3.027  -1.191  -6.987                       H
HETATM   13  H           0      -2.127   1.603  -5.928                       H
HETATM   14  H           0      -2.334   0.779  -3.108                       H
HETATM   15  H           0      -3.817   0.247  -2.352                       H
END
CONECT    1    3    2   14   15
CONECT    2    1
CONECT    3    1    4    5
CONECT    4    3
CONECT    5    3    6    9
CONECT    6    5    7    8   12
CONECT    7    6
CONECT    8    6
CONECT    9    5   10   11   13
CONECT   10    9
CONECT   11    9
CONECT   12    6
CONECT   13    9
CONECT   14    1
CONECT   15    1

mol2.pdb
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33	TITLE Molecule Name REMARK 1 File created by GaussView 5.0.9 HETATM 1 H 0 -2.409 -0.788 -2.338 H HETATM 2 C 0 -2.958 -0.067 -2.907 C HETATM 3 H 0 -3.027 -1.191 -6.987 H HETATM 4 O 0 -3.739 -1.861 -4.248 O HETATM 5 C 0 -3.408 -0.697 -4.233 C HETATM 6 H 0 -2.127 1.603 -5.928 H HETATM 7 N 0 -3.427 0.073 -5.380 N HETATM 8 C 0 -3.804 -0.532 -6.661 C HETATM 9 H 0 -4.713 -1.084 -6.541 H HETATM 10 H 0 -3.948 0.238 -7.390 H HETATM 11 C 0 -3.090 1.494 -5.474 C HETATM 12 H 0 -3.824 1.996 -6.069 H HETATM 13 H 0 -3.073 1.922 -4.493 H HETATM 14 H 0 -2.334 0.779 -3.108 H HETATM 15 H 0 -3.817 0.247 -2.352 H END CONECT 1 2 CONECT 2 1 5 14 15 CONECT 3 8 CONECT 4 5 CONECT 5 2 4 7 CONECT 6 11 CONECT 7 5 8 11 CONECT 8 3 7 9 10 CONECT 9 8 CONECT 10 8 CONECT 11 6 7 12 13 CONECT 12 11 CONECT 13 11 CONECT 14 2 CONECT 15 2

mol2.pdb

TITLE      Molecule Name
REMARK   1 File created by GaussView 5.0.9
HETATM    1  H           0      -2.409  -0.788  -2.338                       H
HETATM    2  C           0      -2.958  -0.067  -2.907                       C
HETATM    3  H           0      -3.027  -1.191  -6.987                       H
HETATM    4  O           0      -3.739  -1.861  -4.248                       O
HETATM    5  C           0      -3.408  -0.697  -4.233                       C
HETATM    6  H           0      -2.127   1.603  -5.928                       H
HETATM    7  N           0      -3.427   0.073  -5.380                       N
HETATM    8  C           0      -3.804  -0.532  -6.661                       C
HETATM    9  H           0      -4.713  -1.084  -6.541                       H
HETATM   10  H           0      -3.948   0.238  -7.390                       H
HETATM   11  C           0      -3.090   1.494  -5.474                       C
HETATM   12  H           0      -3.824   1.996  -6.069                       H
HETATM   13  H           0      -3.073   1.922  -4.493                       H
HETATM   14  H           0      -2.334   0.779  -3.108                       H
HETATM   15  H           0      -3.817   0.247  -2.352                       H
END
CONECT    1    2
CONECT    2    1    5   14   15
CONECT    3    8
CONECT    4    5
CONECT    5    2    4    7
CONECT    6   11
CONECT    7    5    8   11
CONECT    8    3    7    9   10
CONECT    9    8
CONECT   10    8
CONECT   11    6    7   12   13
CONECT   12   11
CONECT   13   11
CONECT   14    2
CONECT   15    2

AMBER建模工具所给构型与其他来源构型的匹配

纤维素二聚体: amber结构与晶体结构
脂分子: amber结构与下载结构

有待改进

或许无须手动指定锚点原子, 而是自动查找测试所有可能的锚点原子, 这样就可以自动处理原子数很多的情况. 具体处理速度待考察.
匹配时仅仅考虑原子间的连接, 没有考虑具体坐标和构型, 无法区分不同的构象.
匹配分两步进行, 首先匹配非氢原子, 然后将氢原子按出现顺序分配到所连重原子上, 因此无法区分不同异构体的氢, 如顺反式.
上面两个缺陷或许可以使用骨架匹配后, 进行最小二乘叠合, 查找最小RMSD的方法解决.