R-group Replacement 实战：如何系统替换取代基

先导优化的日常就是在保持核心不变的前提下，系统地探索取代基空间。这个过程可以做得很随意（想到什么试什么），也可以做得很系统——后者的效率高得多。

第一步：R 基团分解，看清现状

from rdkit import Chem
from rdkit.Chem import rdRGroupDecomposition as rgd
import pandas as pd

core = Chem.MolFromSmarts("c1ccc2c(c1)[nH]c(n2)[*:1]")   # 带连接点的核心
mols = [Chem.MolFromSmiles(s) for s in df["smiles"]]
mols = [m for m in mols if m is not None]

res, unmatched = rgd.RGroupDecompose([core], mols, asSmiles=True, asRows=False)
rg_df = pd.DataFrame(res)
rg_df["pIC50"] = df["pIC50"].values[:len(rg_df)]

print(f"分解成功 {len(rg_df)}，未匹配 {len(unmatched)}")
print(rg_df.head())

# 看每个位点已经试过哪些取代基
for col in rg_df.columns:
    if col.startswith("R"):
        print(f"\n{col}: 已试过 {rg_df[col].nunique()} 种")
        # 按活性排序，看哪些取代基表现好
        top = rg_df.groupby(col)["pIC50"].agg(["mean", "count"])
        print(top.sort_values("mean", ascending=False).head(5))

这张表本身就很有信息量：它直接显示每个位点试过什么、哪些有效、哪些位点探索得还不够。

第二步：系统枚举候选取代基

from rdkit.Chem import AllChem

# 常用取代基库（按性质分组，便于系统探索）
R_GROUPS = {
    "small_alkyl":   ["[*:1]C", "[*:1]CC", "[*:1]C(C)C", "[*:1]C1CC1"],
    "halogen":       ["[*:1]F", "[*:1]Cl", "[*:1]Br", "[*:1]C(F)(F)F"],
    "polar_donor":   ["[*:1]O", "[*:1]N", "[*:1]NC", "[*:1]C(=O)N"],
    "polar_acceptor":["[*:1]OC", "[*:1]C#N", "[*:1]S(=O)(=O)C"],
    "acidic":        ["[*:1]C(=O)O", "[*:1]c1nnn[nH]1", "[*:1]S(=O)(=O)N"],
    "basic":         ["[*:1]N1CCOCC1", "[*:1]N1CCNCC1", "[*:1]CCN(C)C"],
    "aromatic":      ["[*:1]c1ccccc1", "[*:1]c1ccncc1", "[*:1]c1ccc(F)cc1"],
}

def enumerate_analogs(core_smiles, position, r_list):
    """把 R 基团接到核心的指定位点"""
    out = []
    for r in r_list:
        # 用反应模板组合（实践中需处理连接点编号）
        smi = core_smiles.replace(f"[*:{position}]", r.replace("[*:1]", ""))
        m = Chem.MolFromSmiles(smi)
        if m is not None:
            out.append(Chem.MolToSmiles(m))
    return out

all_analogs = []
for group_name, r_list in R_GROUPS.items():
    analogs = enumerate_analogs(core_smiles, 1, r_list)
    for a in analogs:
        all_analogs.append({"smiles": a, "r_group_type": group_name})

analog_df = pd.DataFrame(all_analogs).drop_duplicates("smiles")
print(f"枚举出 {len(analog_df)} 个类似物")

第三步：MMP 分析——从数据学习有效替换

# mmpdb 从真实数据中挖掘「哪些替换实际改善了什么性质」
pip install mmpdb

mmpdb fragment compounds.smi -o frag.fragdb --num-jobs 8
mmpdb index frag.fragdb -o mmp.mmpdb
mmpdb loadprops -p properties.csv mmp.mmpdb

# 查询某个替换的统计效果
mmpdb transform --smiles "目标分子" --property pIC50 --property logD mmp.mmpdb

# 用 Python 分析 MMP 结果
import sqlite3
import pandas as pd

conn = sqlite3.connect("mmp.mmpdb")
query = """
SELECT rule.smiles1, rule.smiles2,
       AVG(pair.delta) AS avg_delta,
       COUNT(*) AS n_pairs,
       STDEV(pair.delta) AS std_delta
FROM ... -- 具体表结构参考 mmpdb 文档
GROUP BY rule.id HAVING n_pairs >= 5
ORDER BY avg_delta DESC
"""

# 关键洞察：
#   n_pairs 大且 std_delta 小的规则 → 效果稳定，值得优先试
#   n_pairs 大但 std_delta 大 → 效果依赖上下文，不可靠
#   n_pairs 小 → 统计不足，仅供参考

「效果稳定性」比「平均效果」更重要：一个平均提升 0.5 log 但标准差 1.2 的替换，实际上是不可预测的；而平均提升 0.3 log、标准差 0.2 的替换更值得信赖。

第四步：多参数优先级排序

import numpy as np
from rdkit.Chem import Descriptors, Crippen, QED

def desirability(v, low, high, direction="range"):
    if direction == "range":
        c, w = (low + high) / 2, max((high - low) / 2, 1e-9)
        return float(np.exp(-((v - c) / w) ** 2))
    if direction == "lower":
        return float(1 / (1 + np.exp((v - high) / (0.1 * abs(high) + 1e-9))))
    return float(1 / (1 + np.exp(-(v - low) / (0.1 * abs(low) + 1e-9))))

def score_analog(smi, pred_activity, pred_activity_std):
    m = Chem.MolFromSmiles(smi)
    if m is None:
        return 0.0
    scores = {
        "activity": desirability(pred_activity, 7.0, 10.0, "higher"),
        "logp":     desirability(Crippen.MolLogP(m), 1.5, 3.5, "range"),
        "tpsa":     desirability(Descriptors.TPSA(m), 50, 110, "range"),
        "mw":       desirability(Descriptors.MolWt(m), 300, 480, "range"),
        "qed":      QED.qed(m),
    }
    weights = {"activity": 2.0, "logp": 1.0, "tpsa": 1.0, "mw": 0.8, "qed": 0.5}
    # 几何平均：任一项差则整体差
    prod = np.prod([scores[k] ** weights[k] for k in scores])
    base = float(prod ** (1 / sum(weights.values())))
    # 不确定性惩罚
    return base * float(np.exp(-0.5 * pred_activity_std))

analog_df["mpo_score"] = [
    score_analog(s, a, sd)
    for s, a, sd in zip(analog_df["smiles"],
                        analog_df["pred_pIC50"],
                        analog_df["pred_std"])
]
top = analog_df.nlargest(50, "mpo_score")

第五步：多样性挑选与实验设计

# 不要只挑分数最高的 —— 要覆盖不同的取代基类型
from rdkit.SimDivFilters import MaxMinPicker
from rdkit.Chem import rdFingerprintGenerator

gen = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=2048)
fps = [gen.GetFingerprint(Chem.MolFromSmiles(s)) for s in top["smiles"]]

picker = MaxMinPicker()
picked_idx = picker.LazyBitVectorPick(fps, len(fps), 20, seed=42)
selected = top.iloc[list(picked_idx)]

# 检查覆盖面：各类取代基都有代表吗？
print(selected["r_group_type"].value_counts())

# 实验设计原则：
#   1) 一次只改一个位点，便于解读 SAR
#   2) 覆盖不同的性质方向（大小、极性、酸碱性）
#   3) 包含 1~2 个「预期会失活」的对照，验证 SAR 假设
#   4) 包含已知活性分子作为实验内标

第 3 条常被忽略但很重要：如果你的 SAR 假设是「这个位点需要氢键供体」，那就该做一个把供体去掉的分子来验证。全做「预期有效」的分子，学不到机制信息。

延伸资源

• 概念：042《R-group Decomposition》；工具：281《Cresset Spark》、291《Boltz Web Demo》；
• 上一步：352《Scaffold Hopping 实战》；MPO 决策见「决策与监管」模块。