APL中的扩展连通性指纹

# APL 中的扩展连通性指纹 来源：https://butwhyisthat.substack.com/p/ecfp 引用 *Rogers 和 Hahn (https://doi.org/10.1021/ci100050t)* 论文 **1** 是每个化学信息学从业者的必经之路。本文试图忠实地实现扩展连通性指纹（ECFP）。我不会费力去解释 ECFP 为何如此流行以及其具体原理。每年大概有十几篇论文中，ECFP 的基线表现似乎好得可疑 **。本文有两个目的：1. 粗略解释 ECFP 的工作原理；2. 在 APL（以及 NumPy 中验证其正确性）中实现它。我相信大多数人要么喜欢 APL，要么讨厌它，要么完全不了解它。所以，我将把这门语言的**学习 (https://www.aplwiki.com/)**部分 (https://www.dyalog.com/getting-started.htm)*留给读者作为*练习 (https://tryapl.org/)*，嘿嘿 (https://tutorial.dyalog.com/)。我阅读并学习了*《精通 Dyalog APL (https://www.dyalog.com/uploads/documents/MasteringDyalogAPL.pdf)》*（我认为这本书写得很好且易于理解）中的内容。在我继续啰嗦之前，我想向读者推荐这些优秀的资源，它们已经很好地解释了 ECFP： 1. https://depth-first.com/articles/2019/01/11/extended-connectivity-fingerprints/ 2. https://cdpkit.org/cdpl_python_cookbook/descr/fingerprints/ecfp.html 3. https://www.blopig.com/blog/tag/extended-connectivity-fingerprints/ 4. https://medium.com/@musicalchemist/from-theory-to-code-a-deep-dive-into-molecular-extended-connectivity-fingerprints-ecfps-with-da1ed436925e 5. https://chemicbook.com/2021/03/25/a-beginners-guide-for-understanding-extended-connectivity-fingerprints.html 6. https://docs.chemaxon.com/latest/fingerprints_extended-connectivity-fingerprint-ecfp.html ECFP 是一种**拓扑指纹 (https://www.blopig.com/blog/2022/06/exploring-topological-fingerprints-in-rdkit/)** **2**。论文开头提到： > [...] ECFP 是通过 Morgan 算法的一种变体衍生出来的 [...] 那么，什么是 **Morgan 算法？** **3**？用通俗的话说，这是一种旨在为化学文摘社（CAS (https://www.cas.org/)）开发基于计算机的登记系统的方法，用于唯一标识化学结构（*在* *20 世纪 60 年代*！！！）。该算法试图为每个二维化学结构分配一个唯一的字母数字描述。  所以，在 60 年代中期（我说得好像我那时在场一样），在披头士狂热、民权运动、暗杀、冷战之间，而且肯定是在 SMILES *之前*，存储化学结构最常见的格式是连接表。连接表是一个有序列表，作为化学结构二维投影的明确的机器描述。它通过描述每个非氢原子的原子符号（或“值”）以及它与图中其他原子的确切键合连接来定义结构。为了捕捉这些结构关系，连接表被组织成五个特定的列表：*FROM ATTACHMENT*（记录编号最小的连接原子）、*RING CLOSURE*（记录形成环的键）、*NODE VALUE*（记录原子符号）、*LINE VALUE* 列表（编码键类型）和 *MODIFICATIONS* 列表（描述特殊属性，如离子电荷、同位素质量或异常化合价）。 敏锐的读者可能已经看出这种表示法在搜索时可能存在的问题，但如果你还不清楚，让我们考虑这个抽象分子的一个任意例子：  具有 5 个原子和 4 个键的任意分子图。你认为应该从 A、B、C、D 还是 E 开始制表？以及这些原子的顺序应该如何？考虑如果我们天真地尝试两种不同的编号方法会发生什么。下面是将它们的 *FROM* 列表并排放置的结果：  同一个图的两种编号方法给出了完全不同的连接表 **4**。你应该存储哪一个？如何选择要存储的？现在，我要说明的是，我只展示了两种图排序，但严格来说，排序空间看起来更像这样 **5**：  希望这能说服你：只为了存储和搜索一个分子，或者为了枚举所有可能的排序并检查是否有任何一个与查询分子匹配，而存储所有这些图是非常不方便的。无论哪种方式，效率都极其低下。如果你仍然不信，让我们看看基于节点数的图的不同枚举数量。这个函数是 *阶乘*（!）  用数字来表明这种情况有多荒谬 **6**：存储*仅一个*具有 27 个原子的*分子*的所有排列，所需的数据存储量将比当前整个数字世界 **7** 的存储量大约多出 50,000 倍。所以这种暴力枚举搜索的想法至少可以说是愚蠢的。现在，*如果*有一种方法可以接受任意连接表/图，并始终以相同的方式为给定分子排列或排序它，那么问题就解决了。这正是 Morgan 所做的。我现在将停止使用表的类比，而是将其视为一个图论问题（仅仅是因为它听起来 *时髦* 和 *酷*）。这个问题可以被视为 *图规范化* 问题（即，给定一个图，输出其顶点的规范标记，该标记仅取决于图的结构）或 *图同构* 问题（即，给定两个图，判断它们在重新标记后是否相同）。规范化是两者中“*更强 (https://math.stackexchange.com/questions/2543452/what-are-the-strong-and-weak-in-mathematics)*”的一个：如果你能规范化，你就可以通过规范化两个图并比较它们来免费判断同构性。无论如何，对于一个给定的图，存在一个能规范地对其节点进行排序的函数，这是解决我们问题的一个机会。 但在我们对此过于兴奋之前，让我们实际看看 Morgan 是如何提出这样做的。 Morgan 算法简单地迭代执行以下步骤： 1. 为每个重原子赋予一个初始连通性值，称为***EC 值***，代表***扩展连通性**，*是的，这正是 ECFP 中的前两个词。这一点很重要：我们创建一个变量***k***，它是分子中唯一 EC 值的数量。例如，如果 EC 值列表是`[1,2,3,4,4]`，那么***k***是`4`。对于懂 Python 的人来说，这很简单： 2. 对于每个原子，计算其*更新后*的 EC 值，即其邻居 EC 值的总和。将不同 *更新后* 值的数量称为***k’***。 3. 如果***k’*****\>*****k***，用*更新后*的值替换*旧* EC 值，将***k***设置为***k’***，并重复步骤 2。 4. 如果***k’*****\≤****k**，退出循环，即停止。 Morgan 实际使用的恒定值是前一次迭代的值，即上次***k***严格增加时的值。当过程终止时，终止迭代的值会被丢弃（*微妙的伏笔，嘿嘿*）。你保留的是上次增加***k***的那组值。Morgan 称之为“分配给节点的最后一组值”。这就是所有下游的破平局（部分排序、规范编号、连接表本身）所依赖的基础。 *困惑吗？* 好吧，看看咖啡因上的迭代过程：  迭代 0 只是重原子的度数。迭代 1 传递/将结构上下文向外扩散：现在原子的值取决于其邻居的连接程度。到迭代 3 时，在“不同距离处看起来不同”的原子（例如靠近氧原子、靠近两个氮原子、在稠环内部等）具有非常不同的值。当进一步细化不再提供新的等价类时，过程停止，只返回前一次迭代的值。好的一面是，只要图中两个原子之间存在真正的结构差异，Morgan 通常会将其揭示出来。例如，萘清晰地划分成三个等价类：桥头碳、α-碳、β-碳，Morgan 在一次迭代后就找到了所有三个类：  按 Morgan 类着色的萘原子 我们可以对我们的简单抽象分子做同样的事情，以使其更清楚。  我们可以让迭代继续，但很明显，在这种情况下，唯一值的数量不会改变。在你继续阅读之前，请猜猜看：***k*** 随着迭代次数的增加而不变是一个普遍趋势吗？  如上所示，抽象图 AB(E)CD 的节点在 Morgan 算法 0 到 25 次迭代中的 EC 值。如预期，唯一值 ***k*** 在迭代 0 后没有变化。 好吧，如果你猜是*假的*，那么你完全正确！事实上，我现在将展示一个非常简单的反例，以证明结构不同的节点可能会*意外地*求和得到完全相同的数字，从而使停止条件*不充分*。  k 的暂时暂停或下降可能仅仅是一次数值碰撞。这并不能从数学上保证算法已经完成了不同类的寻找，因为 k 绝对有可能在后续迭代中再次增加。颜色代表唯一值。 这似乎暗示了一个更普遍的局限性。Morgan 的更新规则通过简单地将邻居值的多重集（一种花哨的说法，即该原子的所有已知/考虑的信息）加起来，将其简化为一个数字。而求和是一个多重集的有损哈希（一种花哨的说法，即两个不同的数字对可以有相同的和）。因此，存在一些原子对在结构上是可区分的（它们的邻居多重集不同），但 Morgan 会给它们分配相同的值。在图论文献中，这一类迭代的颜色细化方案被称为 *Weisfeiler–Leman 算法*（1-WL）**8**。带求和的 Morgan 算法严格*弱于* 1-WL，而 1-WL 本身也有盲点。下面是一个经典例子：两个六顶点的 2-正则图，1-WL 无法区分它们：  1-WL 无法区分的两个六节点 2-正则图 两个图中的每个节点度数均为 2。两个图中每个节点的邻居都形成多重集`{2,2}`。颜色直方图在每次迭代中都一致。但这两个图显然是不同的。左边的图是不连通的。1-WL 看不到这种差异，Morgan 也看不到！ 将 Morgan 算法应用于高度对称的结构也很有趣。例如，这里 ***k*** 保持为 1。例如，环己烷的六个碳是对称的，即没有一个规范过程能仅基于结构将其中一个称为“位置 1”。这种对称性必须在编号阶段被任意打破（在原论文中，通过选择最低的输入索引），在这一点上基本上是随机的。  在完全对称的系统（正则图）上，Morgan 的算法完全不起作用。这显然是预期的，因为无论你如何对这些节点排序，表示总是相同的。平心而论，对于 Morgan 为其构建的 CAS 注册系统来说，这些实际上都不太重要。真实的药物类分子不对称性足够强，以至于基于和的恒定值几乎总能将它们分开，而 Morgan 通过一堆其他破平局方法（节点值 = 原子符号，线值 = 键类型，对端原子的前瞻，在切割键上进行子图分割）来固定剩余部分。 我希望现在你对 Morgan 算法的作用和工作方式有了一些直觉。我可以展示并详细讨论许多边缘情况和有趣的结构，但让我们试着回到本文实际讨论的内容：ECFP。 现在，如果你带着指纹的角度重新阅读最后一节，有趣的一点不是终止准则。而是 Morgan 算法*沿途*生成的东西。每次迭代为每个原子产生一个整数，这个整数总结了半径***r*** 内围绕该原子的所有信息。每个半径每个原子一个整数。如果你贪心一点，你会把这些整数打包起来，并将其用作描述分子的一个*特征包*！这本质上就是 Rogers 和 Hahn 所做的。来自 ECFP 论文： > “ECFP 算法对标准 Morgan 算法做了几处修改。首先，ECFP 生成在预定的迭代次数后终止，而不是在达到标识符唯一性后终止。初始原子标识符以及每次迭代后的所有标识符都被收集到一个集合中；正是这个集合定义了扩展连通性指纹。ECFP 算法不是丢弃中间原子标识符，而是保留它们。” 实际上不多，只有三个修改： 1. **预先固定迭代次数。** 对于 ECFP_*N*，我们执行 *N/2* 次迭代，其中 *N* 是捕获的最大片段（以键计）的直径。每次迭代增加一个键的半径，直径是半径的两倍。三次迭代得到 ECFP_*6*；两次迭代得到 ECFP_*4*；以此类推。 2. **保留所有中间结果。** 每次迭代中计算的每个原子不变量都放入一个列表中。这个列表，去重