AI Drug MapAI 药物发现地图
350

在分子生成中加入 ADMET 约束

在分子生成中加入 ADMET 约束，让产出的分子从一开始就考虑成药性。这篇给出接入自训 QSAR 模型的完整方法与不确定性护栏的设计。

预计 12 分钟 更新于 2026.07.31 实战流程

只优化活性的生成模型会产出「结合很强但完全不能成药」的分子。把 ADMET 预测接进奖励函数，让模型在生成时就考虑溶解度、通透性、代谢稳定性、hERG 风险——这是让生成结果真正可用的关键一步，也是最容易出错的一步。

接入自训 QSAR 模型

# REINVENT4 支持外部模型作为打分组件

[[stage.scoring.component]]
[stage.scoring.component.ExternalProcess]
[[stage.scoring.component.ExternalProcess.endpoint]]
name = "hERG_risk"
weight = 1.0
params.executable = "/path/to/python"
params.args = "/path/to/admet_scorer.py --model herg"
transform.type = "reverse_sigmoid"    # 风险越低越好
transform.low = 4.5                    # pIC50 < 4.5 视为安全
transform.high = 6.0                   # pIC50 > 6.0 视为高风险
transform.k = 0.5
#!/usr/bin/env python3
"""admet_scorer.py —— REINVENT4 外部打分组件
读 stdin 的 JSON {"smiles": [...]}，输出 {"version":1,"payload":{"predictions":[...]}}
"""
import sys, json, argparse
import numpy as np
import joblib
from rdkit import Chem
from rdkit.Chem import rdFingerprintGenerator, DataStructs, Descriptors

gen = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=2048)
DESC = [Descriptors.MolWt, Descriptors.MolLogP, Descriptors.TPSA,
        Descriptors.NumHDonors, Descriptors.NumHAcceptors,
        Descriptors.NumRotatableBonds]

def featurize(smiles):
    X = []
    for s in smiles:
        m = Chem.MolFromSmiles(s)
        if m is None:
            X.append(np.zeros(2048 + len(DESC)))
            continue
        arr = np.zeros(2048, dtype=np.uint8)
        DataStructs.ConvertToNumpyArray(gen.GetFingerprint(m), arr)
        X.append(np.concatenate([arr, [f(m) for f in DESC]]))
    return np.nan_to_num(np.array(X))

def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--model", required=True)
    a = ap.parse_args()

    data = json.load(sys.stdin)
    smiles = data["smiles"]

    # 加载集成模型（用于估计不确定性）
    models = joblib.load(f"models/{a.model}_ensemble.pkl")
    X = featurize(smiles)
    preds = np.stack([m.predict(X) for m in models])
    mean, std = preds.mean(axis=0), preds.std(axis=0)

    # 关键：不确定性高时把预测拉向「保守」方向
    # hERG 是风险指标，不确定时应假设风险较高
    conservative = mean + 1.0 * std

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

if __name__ == "__main__":
    main()

核心问题：模型会钻模型的空子

这是接入 QSAR 组件最危险的地方。生成模型的目标是最大化分数，它会主动搜索预测模型的盲区——那些落在训练分布之外、模型给出虚高分数的区域。结果是生成一批「预测很好但实际完全未知」的分子。

# 应对一：不确定性惩罚（最重要）
def uncertainty_aware_score(mean, std, threshold=6.0, penalty_weight=1.5):
    """预测好但不确定性高的分子，分数要打折"""
    base = 1.0 / (1.0 + np.exp((mean - threshold)))   # 越低越好的性质
    confidence = np.exp(-penalty_weight * std)        # 不确定性越大越接近 0
    return base * confidence

# 应对二：适用域约束 —— 限制与训练集的距离
from rdkit import DataStructs

def applicability_domain_score(smiles, train_fps, gen, k=5, min_sim=0.35):
    """与训练集最近 k 个分子的平均相似度"""
    scores = []
    for s in smiles:
        m = Chem.MolFromSmiles(s)
        if m is None:
            scores.append(0.0); continue
        fp = gen.GetFingerprint(m)
        sims = DataStructs.BulkTanimotoSimilarity(fp, train_fps)
        topk = sorted(sims, reverse=True)[:k]
        avg = float(np.mean(topk))
        # 低于阈值直接归零 —— 明确告诉模型「这里我不知道」
        scores.append(avg if avg >= min_sim else 0.0)
    return np.array(scores)

配置多个 ADMET 组件

[stage.scoring]
type = "geometric_mean"      # 任一项差则总分差

# 溶解度（越高越好）
[[stage.scoring.component]]
[stage.scoring.component.ExternalProcess]
[[stage.scoring.component.ExternalProcess.endpoint]]
name = "logS"
weight = 0.8
params.args = "/path/to/admet_scorer.py --model solubility"
transform.type = "sigmoid"
transform.low = -5.0
transform.high = -3.0
transform.k = 0.5

# 微粒体清除率（越低越好）
[[stage.scoring.component]]
[stage.scoring.component.ExternalProcess]
[[stage.scoring.component.ExternalProcess.endpoint]]
name = "HLM_CLint"
weight = 0.8
params.args = "/path/to/admet_scorer.py --model hlm"
transform.type = "reverse_sigmoid"
transform.low = 10.0
transform.high = 50.0
transform.k = 0.4

# 适用域（关键护栏）
[[stage.scoring.component]]
[stage.scoring.component.ExternalProcess]
[[stage.scoring.component.ExternalProcess.endpoint]]
name = "in_domain"
weight = 1.5                 # 权重要高
params.args = "/path/to/domain_scorer.py"

用简单规则做兜底

在模型不可靠的区域，基于物理化学的规则比模型更稳健：

# 这些规则不依赖训练数据，不会被钻空子
[[stage.scoring.component]]
[stage.scoring.component.TPSA]
[[stage.scoring.component.TPSA.endpoint]]
name = "TPSA"
weight = 0.6
transform.type = "double_sigmoid"
transform.low = 40.0
transform.high = 130.0        # 口服吸收的经验上限
transform.coef_div = 130.0

[[stage.scoring.component]]
[stage.scoring.component.SlogP]
[[stage.scoring.component.SlogP.endpoint]]
name = "cLogP"
weight = 0.8
transform.type = "double_sigmoid"
transform.low = 1.0
transform.high = 4.0          # 过高关联 hERG、代谢不稳定、溶解度差
transform.coef_div = 4.0

[[stage.scoring.component]]
[stage.scoring.component.NumHBD]
[[stage.scoring.component.NumHBD.endpoint]]
name = "HBD"
weight = 0.5
transform.type = "reverse_sigmoid"
transform.low = 2.0
transform.high = 5.0

cLogP 是最有价值的单一约束之一：它与 hERG 风险、代谢不稳定性、溶解度差、非特异性结合都相关。控制住 logP，很多下游问题会同时改善。

生成后的验证

import pandas as pd
import numpy as np

gen_df = pd.read_csv("stage3_1.csv")

# 1) 检查各 ADMET 组件的分布是否真的改善了
for col in [c for c in gen_df.columns if c.startswith("raw_")]:
    print(f"{col}: {gen_df[col].mean():.3f} ± {gen_df[col].std():.3f}")

# 2) 关键检查：高分分子是否落在适用域内
high_score = gen_df.nlargest(100, "Score")
print(f"top100 的适用域分数均值: {high_score['raw_in_domain'].mean():.3f}")
# 如果很低，说明模型正在钻空子 —— 需要加大适用域权重

# 3) 与已知活性分子的相似性分布
# 全都很不相似 → 可能偏离了有效化学空间

# 4) 人工审查（不可替代，见 217）
import mols2grid
mols2grid.display(high_score, smiles_col="SMILES",
                  subset=["img", "Score", "raw_hERG_risk", "raw_SA"],
                  n_cols=5)

常见坑与提示

  • 生成模型会主动搜索 QSAR 模型的盲区，必须加不确定性惩罚与适用域约束；
  • 不确定时把预测拉向保守方向（风险指标假设更高）；
  • 基于物理化学的规则（cLogP、TPSA、HBD）不会被钻空子，是可靠兜底；
  • 生成后必须检查高分分子是否落在适用域内。

延伸资源