子结构搜索在化学信息学中非常常见，利用这个算法，可以从数据库中检索到感兴趣的化合物、可以用来构建基于片段的 QSAR 模型、用来检查化合物中是否存在特定（警示）子结构、甚至用来匹配化学反应规则等等。由于我们可以将化合物分子看成是一个由点和边构成的图（Graph），子结构搜索就成了子图识别的过程。因此我们需要使用一些子图识别的算法来解决子结构搜索的问题，我们一般把子结构认为问问询对象（query），把所搜索的化合物看做目标对象（target）。而其本质又则是图同构（graph isomorphism）问题，即你怎么判断两个图的拓扑结构完全一模一样。

NP 问题

子结构搜索是个 NP（“多项式非确定性”）问题，也就是随着图变得复杂，计算量也会猛增，若没有好的算法进行优化会非常可怕。试想要将问询对象匹配到目标对象，无非就是寻找所有映射对象 M（长度等于问询对象长度 n1），使得问询对象图与映射后的对象（子图）的拓扑结构完全一致。那么理论上就有排列组合 P(n1, n2) 种可能。

比如从一个20个原子的化合物中找6个原子的子结构则可以有27907200种排列（匹配）方式。可以这样认为，假设 query 化合物的六个原子分别是 A1 = [0,1,2,3,4,5]，那匹配过程就是从 target 分子（A2）中挑选 6 个原子，与 A1 中的原子一一匹配，那么 M 的理论数量自然就是全排列情况。当原子数量小时，我们甚至可以使用枚举的方法，比如从[N:0][C:1]=[O:2]中搜索C=O，那么枚举所有两个原子的排列：0,1; 0,2; 1,2; 1,0; 2,0; 2,1共六种可能，不难发现C=O对应1,2。

我们固然可以潜意识的知道可以通过一些手段进行限制、即通过剪枝 (pruning) 减少不必要的考虑。比如当我们知道 target 化合物中的 0 号原子（A2[0]）无法与 A1[0]匹配，那么无需再考虑以 0 开头的任何排列情况，因为这一定是无法匹配的。类似的，比如当A2[2,3]，可以与A1[0,1]匹配，但 A2 的其他原子无法再与A1[2]匹配时，也无需再考虑这些。也可以从图的角度思考，可以把子结构匹配理解为图的生长过程，当某个子图状态 s（比如A2[2,3]）可以和对应状态（A1[0,1]）匹配，则下一步考虑如何进一步生长。

一个合理的算法可以让它避免一些，从而大大减少遍历次数与计算时间。

VF2 算法

现在主流的化学子结构搜索方法都是使用 VF2 算法。不过首先先介绍下大多数子结构匹配算法的本质，可以认为是一下几个步骤的重复：

  1. 初始化状态，s0
  2. 从 query 和 target 中各选一个原子（形成原子对）
  3. 判断这对原子是否可以加入匹配集合 M
  4. 将该原子对加入 M，进入下一个状态 s’，执行2
  5. 回到状态 s 执行2（选择下一个原子对）

以下是原论文[1]中的算法伪代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
PROCEDURE Match(s)
    INPUT: an intermediate state s; the initial state s0 has M(s0)=Æ
    OUTPUT: the mappings between the two graphs
    IF M(s) covers all the nodes of G2 THEN
        OUTPUT M(s)
    ELSE
        Compute the set P(s) of the pairs candidate for inclusion in M(s)
        FOREACH (n, m)ÎP(s)
            IF F(s, n, m) THEN
                Compute the state s´ obtained by adding (n, m) to M(s)
                CALL Match(s¢)
            END IF
        END FOREACH
        Restore data structures
    END IF
END PROCEDURE

从上述算法中可知，其实这是个类似于树的搜索，一般是基于深度的搜索，因为无需遍历所有情况，只需找到一个最深（即所有 query 原子都被匹配）即可。

VF2 算法同样也是这几个步骤，只是此过程中运用了一些技巧，可以减少不必要的搜索。

  1. 在选择一个原子对的时候，首先先从“已匹配原子的相邻原子”。 T1、T2（Terminal的意思，1、2 分别对应 query 和 target）集合中挑选。这非常合乎情理，因为我们肯定希望匹配集 M 是一个“生长”的过程，而不是无序随机的选择一对看最终拓扑结构是否一致。
  2. 在判断原子对是否可以加入匹配集合时，除了看“点的属性”，即原子类型时，需要判断 query 中原子的“度”（即相邻原子数量）是否大于 target 中的，如果大于则直接剪枝。很好理解，如果 query 度大于 target，则将来一定无法让 query 上该原子的邻原子都与 target 匹配上。
  3. 在加入该原子对时，同时更新 T1 与 T2，用以下一层（s’）匹配时选择原子对。

我当时有一点点疑惑，增加 T 变量难道不是增加了计算量么？但其实，增加 T 的计算量其实跟“搜索次数”在一个量级，比如原来搜索 10 次是 10 毫秒，现在长到 20 毫秒，搜索 100 次也不过增长到 200 毫秒，即线性或者顶多多项式倍增长，复杂度没有变。但通过这样的计算可以大大减少“不必要的搜索”，可能原来要搜索 100 次，现在只需要 10 次，最终只花了 20 毫秒。

RDKit 代码简介

无论是 OpenBabel、RDKit 还是 CDK 甚至是 ChemAxon 使用的都是 VF2 算法。原来 RDKit 是最模仿 VF2 算法的原始代码的，可能原因是 RDKit 也是用 C++ 写的，所以无需有太多的“理解和翻译”，即可直接“抄袭代码”，不过 2019 年时已经被优化过，因为原 VF2 算法针对有向图，而分子是无向图，因此很多变量与计算都可以省略。

RDKit 的子结构核心代码在 /Code/GraphMol/Substruct/SubstructMatch.cpp 中，可以从 此处 开始读。这里面作者实例了三个类，分别是 AtomLabelFunctor，BondLabelFunctor以及MolMatchFinalCheckFunctor。这三个对象分别是判断函数，用来判断原子、键以及最终的化合物是否完全匹配上。而从 VF2 算法则是在 这行 的vf2_all函数。而具体的 VF2 的算法则在 vf2.hpp。

NetworkX 代码解析

除了可以参考 RDKit 的源码，笔者更推荐研究 NetworkX 中的 vf2 实现方法（点击此处查看 GitHub 中的源码），由于是使用 python 语言加上面向对象的思路，可能会更好理解一些，这里对两个核心函数进行解释。首先要说明代码里的一些术语，G1 对应上述中的 target，G2 对应 query，即“从 G1 中得到一个子图与 G2 完全一样”。core_1 是 G1 上的点到 G2 上的点的映射，比如 {0: 1}，表示 G1 上 0 位的点与 G2 上 1 位的点匹配上，core_2 则是 G2 到 G1 的映射。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 核心的深度优先搜索算法
def match(self):
    # 搜索结束判断：当 core_1 长度等于 G2 的长度时，
    # 即所有 G2 上的点都能从 G1 中（通过 core_1）找到
    if len(self.core_1) == len(self.G2):
        self.mapping = self.core_1.copy()
        yield self.mapping
    else:
        # candidate_paris_iter 即“寻找下一个原子对”
        for G1_node, G2_node in self.candidate_pairs_iter():
            # 判断 G1_node 与 G2_node 即该原子对是否能匹配上，
            # 在子结构搜索中即判断原子是否一致、连接度是否一致等
            if self.syntactic_feasibility(G1_node, G2_node):
                if self.semantic_feasibility(G1_node, G2_node):
                    # Recursive call, adding the feasible state.
                    # 如果一致，则进入下一层递归，因此是基于深度的搜索
                    newstate = self.state.__class__(self, G1_node, G2_node)
                    yield from self.match()
                    newstate.restore()

下面看下如何寻找下一个原子对，术语解释：self.inout_1就是与 G1 中已匹配的（core_1）点相邻的点。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def candidate_pairs_iter(self):
    G1_nodes = self.G1_nodes
    G2_nodes = self.G2_nodes
    min_key = self.G2_node_order.__getitem__

    # 计算 T1 和 T2，即已经匹配上的原子的相邻原子
    T1_inout = [node for node in self.inout_1 if node not in self.core_1]
    T2_inout = [node for node in self.inout_2 if node not in self.core_2]

    # 如果同时都有相邻的原子，G2 选择顺序最小的原子，G1 则遍历所有可能。
    # 这是因为 G2 是 query，因此必然每个原子都需要被匹配上，
    # 因此这里只需要选择一个即可，剩下的原子下次仍然需要匹配
    if T1_inout and T2_inout:
        node_2 = min(T2_inout, key=min_key)
        for node_1 in T1_inout:
            yield node_1, node_2
    else:
        # 当没有相邻原子时，比如开始时，或者 query 分子是“两个分子”时
        # 从 G2 中剩下的原子中选择索引最小的值
        other_node = min(G2_nodes - set(self.core_2), key=min_key)
        # 遍历所有 G1 所有剩下的原子
        for node in self.G1:
            if node not in self.core_1:
                yield node, other_node

SMARTS 匹配

虽然 SMARTS 匹配超出了本文的范围，不过作为化学信息学中非常常见的子结构搜索需求，还是需要说明一下。SMARTS 匹配的核心就在于所需匹配的原子并非明确的真实原子，而是一个“问询原子（QueryAtom）”，比如在SMARTS中可以用“[N,O]”表示一个氮原子或者氧原子的原子。其详细规则可见 Daylight 说明。 SMARTS 匹配与上述 VF2 算法的藕合处就在于 VF2 中的第三步“判断这对原子是否可以加入匹配集合 M”，判断的原子就是“原子匹配”、“键匹配”，对应 RDKit 中的AtomLabelFunctor和BondLabelFunctor两个函数。

[1] L. P. Cordella, P. Foggia, C. Sansone, M. Vent. An Improved Algorithm for Matching Large Graphs