Nat. Commun. | 当前小分子数据集存在显著覆盖偏差！耶拿大学团队新发现

DRUGAI

今天为大家介绍的是来自德国耶拿大学Sebastian Böcker团队的一篇论文。小分子机器学习旨在从分子结构预测化学、生物化学或生物学特性，其应用包括毒性预测、配体结合和药代动力学。最近的一个趋势是开发避免显式领域知识的端到端模型。这些模型假设训练和评估数据中没有覆盖偏差，这意味着数据能代表真实分布。然而，这类模型很少考虑适用性范围。在这项研究中，作者调查了大规模数据集对已知生物分子结构空间的覆盖程度。为此，作者提出了一种基于求解最大公共边子图(MCES)问题的距离度量方法，该方法与化学相似性有很好的一致性。尽管这种方法在计算上很困难，但作者引入了一种结合整数线性规划和启发式边界的高效方法。研究结果表明，许多广泛使用的数据集缺乏对生物分子结构的均匀覆盖，这限制了在这些数据集上训练的模型的预测能力。作者提出了两种额外的方法来评估训练数据集是否偏离已知分子分布，这可能有助于指导未来数据集的创建以提高模型性能。

几十年来，机器学习在生物化学和化学领域已经得到成功应用。作者着眼于从生物相关小分子的分子结构预测其化学、生物化学或生物学特性这一任务。近期的趋势是开发端到端模型，避免通过归纳偏差显式整合领域知识。值得注意的例子包括用于新型抗生素和高毒性小分子的生成模型，以及用于抗生素活性、气味感知和酶-底物预测的分类器。2018年的MoleculeNet论文提供了17个中等到大规模的分子性质预测数据集。这些数据经常用于机器学习，以训练和评估新模型，如图神经网络和graphormers；该论文在5年内获得了2000多次引用。

在化学计量学界，众所周知不应该在模型适用范围之外使用模型。这种情况可以类比于空间偏差，即使用来自某个地理位置的测试(和训练)数据，但声称模型对其他地理位置也同样有效。然而，在训练用于预测分子性质的大规模端到端模型时，这个问题通常被忽视。其他模型虽然在更大的结构数据集上进行预训练，但这无法解决单个分子性质的训练和评估数据中的分布偏差问题。最近有警告指出，机器学习可能导致科学领域的可重复性危机。特别是，MoleculeNet的数据集受到了批评。相比之下，训练一个在评估中表现良好的机器学习模型相对简单，但要得到一个确实有助于解决基本问题的模型则要困难得多。

数据集内部泛化的问题已经得到广泛研究。对于小分子，广泛使用的骨架分割确保了对训练数据中未见过的骨架进行评估。这使作者能够评估模型对新型分子结构的外推能力。然而，这样做并没有考虑到分子性质分布的差异。特别是，小的结构变化可能导致待预测的相关分子性质产生巨大差异，这种现象被称为"活性悬崖"（activity cliff）。最后，通过采用骨架分割得出的所有外推性能结论仍然只适用于数据集本身及其（受限的）化学空间。在仅由脂质组成的数据集上训练和评估的模型可能在不同脂质类别之间具有良好的泛化能力，但没有理由假设它对类黄酮也同样有效。

为了确保模型不会在其适用范围之外使用，必须确保用于训练和评估模型的可用数据是所有感兴趣分子空间的代表性子集。作者认为，如果数据集不够全面，就无法学习所要解决问题的所有方面。然而，即使对于包含数千个样本的数据集，纳入其中的小分子的选择往往远非随机。对于依赖实验测量的数据集，它反而受化合物可获得性（availability）的支配，因此常常受到经济因素的影响。化合物的可获得性取决于化学合成的难度、前体化合物的商业可得性以及合成化学和生物技术的类似考虑因素。化合物的可获得性越低，价格就越高，这种化合物在大规模数据集中出现的可能性就越小。显然，这给训练数据引入了偏差。对于具有实验测量的数据集，某些化合物的不可获得性在近期内不会改变。

要考虑小分子的训练数据分布，需要一些方法来估计分子结构之间的相似性或差异性。不幸的是，这是一个非常复杂的问题，目前有两种方法在处理这个问题，但各有不足：首先，分子指纹允许快速处理大型数据集。然而，基于分子指纹的测量方法已知具有不理想的特征。特别是，测量的距离可能与化学直觉有很大差异。其次，基于计算最大公共(边)子图的方法能更好地捕捉结构相似性的化学直觉，但不幸的是，需要解决计算上的难题。

在这里，作者展示了如何检查分子结构数据集对生物相关小分子结构（为简便起见称为"生物分子结构"）的覆盖情况。作者的方法结合了统一流形近似和投影(UMAP)嵌入和通过最大公共边子图(MCES)计算结构距离。

对生物分子结构分布的探查

作者使用14个分子结构数据库(简称生物分子结构)的组合作为这个空间的代表。这些数据库包含代谢物、药物、毒素和其他具有生物学意义的小分子。据此，这些数据库总共包含718,097个生物分子结构。显然，这个代表性数据集现在和将来都是不完整的；然而，适用性范围的限制可能已经在这个代表性数据集中显现出来。

给定一对分子结构，作者使用它们的最大公共边子图计算距离。为了加快计算速度，作者估计了所有距离的(可证明正确的)下界。作者仅在距离界限最多为所选距离阈值时才进行精确计算，除非另有说明，这个阈值被设定为10。如果下界高于阈值，作者就使用这个界限作为距离估计。同样，如果计算出的精确距离高于阈值，作者就使用阈值代替。作者使用UMAP来在二维图中可视化生物分子结构的全体。

为了避免运行时间过长和图表过于杂乱，作者均匀抽样了20,000个生物分子结构(图1)。在40核处理器上，MCES计算的总运行时间约为15.5天。为了监测抽样的影响，作者从这个集合中又均匀地二次抽样了9次，每次10,000个分子结构。作者在图2中展示了相应的UMAP嵌入，以显示变化。作者观察到抽样确实可能改变UMAP嵌入的整体布局，但整体布局通常惊人地相似。众所周知，这些UMAP嵌入的解释需要谨慎。

图 1

图 2

某些分子结构和化合物类别，特别是某些脂质类别，在UMAP嵌入中形成了离群簇(图1)。为了避免这些分子结构主导UMAP嵌入，作者将它们排除在外。这样剩下了18,096个分子结构，这些结构将在所有后续分析中被考虑。图3展示了生物分子结构的最终图，其中作者根据ClassyFire对化合物类别进行了颜色编码。被排除的分子结构仍然可以在UMAP嵌入中显示。

图 3

上文中，作者使用UMAP来可视化MCES距离。显然，任何其他用于投影以距离形式给出的高维数据的方法(例如，t-SNE、多维缩放或最小生成树)也可以应用。

图 4

图4展示了近视MCES距离的分布。正如预期的那样，大多数距离都很大；然而，对于每个分子结构来说，与另一个分子结构的最小距离通常低于10(图4b)。实际上，作者观察到距离10出现的频率远高于随机情况下的预期。这是由于作者在计算中使用的阈值T=10，特别是由于双重阈值处理：当计算出精确的MCES距离但发现大于T时，作者就用T作为近视距离。鉴于在近视MCES距离的分布中没有这样的凸起(图4a)，作者认为这种阈值处理对计算出的UMAP嵌入的影响很小。

Chari等人表明，从二维UMAP嵌入中推断数据结构时必须高度谨慎。作者强调这对作者正在做的工作并不构成限制：作者已经知道数据是有结构的，即分子结构可以(不)相似，或者属于或不属于相同的化合物类别。作者正在研究的是数据子集在多大程度上是均匀的子样本：如果作者能在二维UMAP嵌入中发现非均匀性，那么在更高维度中它可能也不是均匀子样本。然而，不可否认的是UMAP嵌入远非完美，这意味着图中的小/大距离不一定意味着小/大的MCES距离。UMAP与任何将高维空间投影到平面的方法都有这些限制。显然，UMAP图的布局也有一定的任意性，参见图2。最后，当将新样本插入现有的UMAP嵌入时，这些新样本倾向于被插入到图的现有结构中，而不是生成新的远距离簇或单个点。

最大公共边子图计算

最大公共边子图问题非常适合用来捕捉分子结构之间的高度相似性。两个分子结构之间的MCES距离等于第一个分子图中的边数加上第二个分子图中的边数，再减去MCES中边数的两倍。这个距离与化学直觉很好地吻合：它是从两个图中需要移除的最少边数，使得剩余的图同构（忽略单个节点）。从化学角度来看，可以将其理解为将一个分子转化为另一个分子所需的化学反应数量。分子图是标记图，其中节点与原子类型相关联，边与键级相关联。在比较两个这样的图时，需要考虑这些标记以及如何比较它们，例如双键与芳香键的比较。在计算中不考虑氢原子。

不幸的是，计算MCES是一个可证明的困难问题，人们认为不太可能存在解决该问题的高效算法。为了计算图3中UMAP嵌入的距离，需要快速解决超过1.6亿个NP难的MCES问题实例。作者提出了一种基于计算距离下界的高效实现方法，用于快速识别不相似的结构对。对于可能相似的结构对，作者使用整数线性规划（Integer Linear Programming，ILP）方法计算精确距离。使用ILP计算MCES相比之前基于乘积图中枚举团的方法有两个关键优势：(i) 与基于团的方法相反，当输入相似时ILP往往更快，这恰恰是最有意思的情况；(ii) ILP避免了通过线图建模问题时繁琐的ΔY变换处理。作者的方法是首次将ILP用于化学结构比较。

在实践中，作者为距离阈值（比如T=10个边修改）判断边界是否能保证真实距离至少为T；在这些情况下，作者使用边界作为真实距离的近似值。只有当没有边界能保证距离至少为T时，作者才执行精确算法并报告结果，这时将T用作上界。作者将得到的距离称为近视MCES距离。这种两步法有两个优点：(a) 当MCES距离较小时，ILP通常较快，而对于较大距离，运行时间可能会非常长。因此作者的两步方法明确有效地排除了大多数耗时的精确计算。(b) 考虑一对分子结构时，知道它们之间的MCES距离是2还是8显然很有意义。但作者认为，知道MCES距离是42还是48大多是无关紧要的：在这两种情况下，两个分子结构都是高度不相似的。这两个数字都可以作为真实距离的合理近似值。特别是，根据UMAP优化问题的表述，较大距离的差异不会导致UMAP嵌入发生实质性变化。

图 5

作者对方法进行了深入评估，使用了20,000个均匀抽样的实例，其中每个实例是一对生物分子结构（图5）。运行时间在一个运行80个并行线程的40核处理器上测量，作者报告了每个线程的运行时间。对于整数线性规划（ILP），20000个实例中有24个在四天的挂钟时间内未能完成。对于这些实例，作者使用计算停止时的时间作为运行时间的代理值。这样，20,000个实例的ILP总运行时间等于234.2天（平均每个实例16.9分钟）。将实例按ILP运行时间排序，作者观察到1.04%的实例占总运行时间的95%以上。

作者首先考虑ILP运行时间与精确MCES距离的关系（图5a，b）。未完成的ILP实例被排除在此图之外。对于75以内的距离，作者观察到距离和运行时间之间存在明显的相关性。例如，29.5%的实例的MCES距离在30以内，但仅贡献了总运行时间的0.91%。对于更大的距离（最高到100），这种相关性变得不那么明显。显然，MCES距离和MCES实例的实际结构都会影响运行时间。对于超过100的距离，结果的解释必须谨慎，因为这些实例较少，对应于少数异常结构。

为了排除通过抽样产生的偏差，作者使用19,994个生物分子结构的所有配对重复了上述分析（参见分子结构抽样部分）。为了避免运行时间激增——在本文的其余部分中也是如此——作者集中在阈值T=10的情况，见图5e，f。完整数据集的总运行时间为1240天，相当于在40核处理器上运行15.5天的挂钟时间。作者观察到基本上与抽样实例相同的模式，但每个实例的平均运行时间比抽样集大5.2倍（536毫秒比102毫秒）。在这里，0.2%的实例占据了总运行时间的一半。在相同阈值下平均运行时间的差异很可能可以归因于ILP求解器设置的不同。作者没有对ILP求解器使用运行时间限制，具体见方法部分。

最后，作者分析了计算界限而不是精确距离时的误差：近视MCES距离可能与精确距离不同；绝对误差是精确距离值与近视MCES距离之间的差值。需要注意的是，近视MCES距离是一个下限，所以绝对误差永远不会是负数。与真实值相比，偏离真实值的情况最为重要。为此，作者研究了相对误差，即将绝对误差除以精确距离值（图5g）。正如预期的那样，如果降低阈值T，更多的距离估计会出现非零误差：具体来说，当T=25时，22.5%的实例误差为零；当T=10时，这个比例降至9.5%。除此之外，作者观察到在10到25之间改变阈值T并不会显著改变相对误差分布的形状。作者发现，对于不同的阈值和精确计算，近视MCES距离的分布（与图4相比）非常相似，只是发生了偏移（图5h）。

最近，RDKit中提供了RASCAL的实现。RASCAL解决了MCES问题，当预定义的相似度阈值被超过时，计算精确的MCES。与作者的实现类似，当不进行精确计算时，计算得到的界限结果可用于近似相似度。与（绝对）近视MCES距离不同，这里使用了Johnson引入的相对相似度度量。为了评估RDKit的RASCAL实现，作者使用了与之前相同的20,000个实例子集。使用默认的相似度阈值0.7，只有280个实例（1.40%）进行了精确计算。值得注意的是，这280个实例中有100个（35.7%）计算失败，原因是超过了产物图的内部大小限制（"maxNumberMatchingBonds"）（98个实例），或者达到了一小时的超时限制（2个实例）。在禁用产物图的内部大小限制后，这98个实例中的68个在不到一小时内完成了精确计算。在29个案例中达到了运行时间阈值，在一个案例中超过了256 GB的可用内存。对于所有279个未超过内存的精确计算实例，每个实例的平均运行时间为7.92分钟，对于超时的实例则以一小时为代表值。

接下来，为了精确计算更多实例，作者将相似度阈值降低到0.5。这样做后，2899个实例（14.5%）被精确计算出来。在这些实例中，423个（14.6%）由于内部大小限制（408个实例）或超时（15个实例）而未能完成计算。当禁用内部大小限制时，408个实例中的186个超过了一小时的时间限制；有5个实例超出了内存限制。在这里，2894个精确计算的实例的平均运行时间为4.47分钟，这不包括因内存使用超限而崩溃的实例。同样，对于超时的实例，使用一小时作为代理时间。最后，通过RDKit的RASCAL计算近似结果比计算MCES边界要快几个数量级，前者的平均运行时间为1.36毫秒。然而，比较RDKit和作者的MCES代码之间的运行时间可能会产生误导，因为这主要是由已编译的C++和解释型Python代码之间的运行时间差异所决定的。

MCES问题存在许多变体,例如寻找最大公共子图(MCS)、连通MCES或该问题的受限变体。比较不同变体结果的质量极具挑战性,因此作者不进行深入评估。Seipp进行了近视MCES距离与使用RDKit的MCS计算和SMARTScompare的运行时间评估,发现近视MCES方法比其中实现的精确MCS方法快几个数量级,有时甚至可以与启发式方法相媲美。这一点值得注意,因为MCES问题通常被认为比MCS问题要困难得多。

对于许多计算任务,需要距离度量是一个度量标准。作者可以轻松证明边界和精确MCES距离都是(伪)度量。对于近视MCES距离,重要的是作者应用双重阈值:一方面,仅在边界低于阈值T时才计算精确距离。另一方面,如果精确距离高于T,作者会将T报告为近视MCES距离。这意味着每当执行精确计算时,报告的距离小于或等于T;而如果只执行边界计算,报告的距离大于或等于T。这样做可以通过约束(8)加快ILP计算。除此之外,它还允许作者证明近视MCES距离确实是一个度量标准。

如果不使用第二个阈值(即,即使精确MCES距离大于T也返回),这不仅会增加ILP的运行时间。此外,由此产生的距离度量不再是一个度量标准。可以使用Floyd-Warshall算法来寻找完整图中的所有最短路径,以强制执行三角不等式。

编译|黄海涛

审稿|王梓旭

参考资料

Kretschmer, F., Seipp, J., Ludwig, M., Klau, G. W., & Böcker, S. (2025). Coverage bias in small molecule machine learning. Nature Communications, 16(1), 554.