在分子生成中加入 Docking Score 约束

把对接分数接进生成的奖励函数，是「基于结构的生成设计」最直接的做法：模型在生成时就考虑分子能否与口袋互补。但对接很慢，且打分函数本身不可靠——这两点决定了实现方式与风险控制。

核心矛盾：速度

REINVENT4 每步生成 128 个分子，跑几百步就是数万次对接。如果每个分子对接需要 10 秒，一轮训练要跑几天。必须做性能优化：

策略	加速	代价
降低 exhaustiveness（8 → 4）	约 2×	姿势质量下降
缩小盒子	1.5～2×	可能切掉正确姿势
多进程并行	核数倍	内存占用
结果缓存	视重复率	无
代理模型	100～1000×	需先积累数据
GPU 对接（GNINA）	数倍	需 GPU

实现：带缓存的对接打分组件

#!/usr/bin/env python3
"""docking_scorer.py —— 带缓存与并行的对接打分"""
import sys, json, os, hashlib, tempfile, subprocess
import numpy as np
from concurrent.futures import ProcessPoolExecutor
from rdkit import Chem
from rdkit.Chem import AllChem

CACHE_FILE = "docking_cache.json"
CENTER = [11.2, 23.4, -5.7]
BOX = [20.0, 20.0, 20.0]

def load_cache():
    return json.load(open(CACHE_FILE)) if os.path.exists(CACHE_FILE) else {}

def canonical_key(smi):
    m = Chem.MolFromSmiles(smi)
    return Chem.MolToSmiles(m) if m else None

def dock_one(smi):
    """返回对接分数（kcal/mol），失败返回 None"""
    mol = Chem.MolFromSmiles(smi)
    if mol is None:
        return None
    try:
        mol = Chem.AddHs(mol)
        p = AllChem.ETKDGv3(); p.randomSeed = 42
        if AllChem.EmbedMolecule(mol, p) != 0:
            return None
        AllChem.MMFFOptimizeMolecule(mol)

        with tempfile.TemporaryDirectory() as td:
            sdf = os.path.join(td, "lig.sdf")
            pdbqt = os.path.join(td, "lig.pdbqt")
            out = os.path.join(td, "out.pdbqt")
            Chem.MolToMolFile(mol, sdf)
            subprocess.run(["mk_prepare_ligand.py", "-i", sdf, "-o", pdbqt],
                           check=True, capture_output=True, timeout=60)
            r = subprocess.run([
                "vina", "--receptor", "receptor.pdbqt", "--ligand", pdbqt,
                "--center_x", str(CENTER[0]), "--center_y", str(CENTER[1]),
                "--center_z", str(CENTER[2]),
                "--size_x", str(BOX[0]), "--size_y", str(BOX[1]),
                "--size_z", str(BOX[2]),
                "--exhaustiveness", "4",       # 生成阶段用低值换速度
                "--num_modes", "1", "--seed", "42",
                "--out", out,
            ], capture_output=True, text=True, timeout=120, check=True)

            for line in r.stdout.split("\n"):
                parts = line.split()
                if len(parts) >= 2 and parts[0] == "1":
                    return float(parts[1])
    except Exception:
        return None
    return None

def main():
    data = json.load(sys.stdin)
    smiles = data["smiles"]
    cache = load_cache()

    todo = [s for s in smiles
            if (k := canonical_key(s)) and k not in cache]
    if todo:
        with ProcessPoolExecutor(max_workers=16) as ex:
            for s, score in zip(todo, ex.map(dock_one, todo)):
                k = canonical_key(s)
                if k:
                    cache[k] = score if score is not None else 0.0
        json.dump(cache, open(CACHE_FILE, "w"))

    scores = []
    for s in smiles:
        k = canonical_key(s)
        v = cache.get(k, 0.0) if k else 0.0
        scores.append(v if v is not None else 0.0)

    json.dump({"version": 1, "payload": {"predictions": scores}}, sys.stdout)

if __name__ == "__main__":
    main()

配置

[[stage.scoring.component]]
[stage.scoring.component.ExternalProcess]
[[stage.scoring.component.ExternalProcess.endpoint]]
name = "docking"
weight = 1.2
params.executable = "/path/to/python"
params.args = "/path/to/docking_scorer.py"
# Vina 分数越负越好，需反向映射
transform.type = "reverse_sigmoid"
transform.low = -11.0        # -11 视为很好
transform.high = -7.0        # -7 视为一般
transform.k = 0.4

更快的方案：代理模型

# 思路：先对接一批分子，训一个快速模型预测对接分数，
#      生成时用代理模型，定期用真实对接校准

import lightgbm as lgb
import numpy as np

class DockingSurrogate:
    def __init__(self, retrain_every=2000):
        self.model = None
        self.X, self.y = [], []
        self.retrain_every = retrain_every
        self.n_since_retrain = 0

    def add_real(self, features, score):
        self.X.append(features); self.y.append(score)
        self.n_since_retrain += 1
        if self.n_since_retrain >= self.retrain_every:
            self.retrain()

    def retrain(self):
        self.model = lgb.LGBMRegressor(n_estimators=500, num_leaves=63,
                                       verbose=-1)
        self.model.fit(np.array(self.X), np.array(self.y))
        self.n_since_retrain = 0

    def predict(self, features):
        if self.model is None:
            return None
        return self.model.predict(features)

# 混合策略：
#   90% 的分子用代理模型（毫秒级）
#   10% 随机抽样做真实对接，用于持续校准与检测代理漂移

对接约束的固有风险

打分函数会被钻空子。这是最严重的问题。Vina 分数偏好大分子和高疏水性分子，生成模型很快会发现这一点，产出一堆又大又油的分子——分数很高，但溶解度极差、选择性差、根本不能成药。

必须同时约束分子大小与疏水性：

# 用配体高效性代替原始分数，或同时加 MW/logP 约束
def ligand_efficiency(score, mol):
    n_heavy = mol.GetNumHeavyAtoms()
    return -score / n_heavy if n_heavy > 0 else 0.0
# LE 对分子大小是归一化的，不会奖励「靠变大刷分」

低 exhaustiveness 的姿势不可靠：为了速度用 exhaustiveness=4，姿势质量下降，分数噪声更大。
构象生成失败会静默返回 0：柔性大的分子可能嵌入失败，被当成「分数为 0」，反而可能被选中。要显式处理失败情况。
对接分数与实测活性相关性弱：即使一切正常，Vina 分数与 IC50 的相关性通常也只有 0.3～0.5。它提供的是方向性引导，不是活性预测。

更稳健的做法

# 用「相互作用指纹匹配度」代替原始对接分数
# 判据不是「分数多高」，而是「是否复现了已知活性分子的关键相互作用」

# 1) 先用共晶配体算出参考相互作用指纹（见 338）
# 2) 生成的分子对接后算指纹
# 3) 用与参考的相似度作为奖励
#
# 优点：不会被「大而油」的分子欺骗，
#      因为它要求形成特定的相互作用模式

常见坑与提示

对接分数偏好大而疏水的分子，必须同时约束 MW 与 logP，或改用配体高效性；
缓存 + 并行 + 低 exhaustiveness 是让对接约束可行的必要优化；
构象生成失败要显式处理，不能静默给 0 分；
用相互作用指纹匹配度代替原始分数，更不容易被钻空子。

延伸资源

生成平台：347《用 REINVENT4 生成新分子》、221《REINVENT4》；对接：336《用 AutoDock Vina 跑 Docking》、337《用 GNINA 做 Rescoring》；
相互作用指纹：338《用 ProLIF 提取相互作用指纹》；其它约束：348《在分子生成中加入 SA Score 约束》、350《在分子生成中加入 ADMET 约束》。