QSAR 入门：从分子结构预测生物活性

QSAR（定量构效关系）用分子结构预测生物活性。它是化学信息学中最古老也最实用的建模任务——现代的「AI 活性预测」本质上就是 QSAR，只是换了更复杂的模型。

建模的基本流程

# 【六个步骤，每一步都可能出错】
#
# 1) 【数据收集与清洗】← 【最关键，也最费时】
# 2) 【分子表示】
# 3) 【数据划分】← 【最容易出错】
# 4) 模型训练
# 5) 【评测】
# 6) 【适用域判断】
#
# 【常见的错误分配】：
#   新手把 80% 时间花在第 4 步（调模型）
#   而实际上第 1、3 步的影响大得多

第一步：数据准备

import pandas as pd
import numpy as np
from rdkit import Chem
from rdkit.Chem.MolStandardize import rdMolStandardize

# ---- 从 ChEMBL 获取数据（见 226）----
# 关键的筛选条件：
#   - 同一个靶点（target_chembl_id）
#   - 【同一种活性类型】（IC50 / Ki / EC50 不要混）
#   - 【同一种测定方法】（如果信息可得）
#   - standard_relation = '='（排除 > 和 <）

def prepare_qsar_data(df):
    # 1) 【转成对数尺度】
    #    活性跨越多个数量级，必须取对数
    df = df[df["standard_units"] == "nM"].copy()
    df["pActivity"] = -np.log10(df["standard_value"] * 1e-9)

    # 2) 【标准化结构】（见 046）
    def standardize(smi):
        m = Chem.MolFromSmiles(smi)
        if m is None:
            return None
        m = rdMolStandardize.Cleanup(m)
        m = rdMolStandardize.FragmentParent(m)     # 去盐
        return Chem.MolToSmiles(m)

    df["std_smiles"] = df["canonical_smiles"].apply(standardize)
    df = df.dropna(subset=["std_smiles"])

    # 3) 【处理重复测量】
    #    同一分子多次测定 → 取中位数（比均值稳健）
    #    【但先检查一致性】：
    grouped = df.groupby("std_smiles")["pActivity"]
    spread = grouped.max() - grouped.min()
    inconsistent = spread[spread > 1.0]     # 差异超过 10 倍
    print(f"【{len(inconsistent)} 个分子的重复测定差异 > 1 个 log 单位】")
    #    → 这些数据可能有问题，考虑剔除

    df = df.groupby("std_smiles", as_index=False).agg(
        pActivity=("pActivity", "median"),
        n_measurements=("pActivity", "count"),
    )

    # 4) 【检查活性分布】
    print(df["pActivity"].describe())
    #    分布过窄 → 建模没有意义
    #    极端值 → 检查是否为异常数据

    return df

# 【删失数据的处理】：
#   "IC50 > 10000 nM" 表示未测到活性
#   → 【不能当作 10000 nM 处理】（会低估）
#   → 也不能直接丢弃（丢失了「不活性」的信息）
#   → 【正确做法】：
#     - 做分类任务时，归为「非活性」类
#     - 做回归时，用删失回归（Tobit 模型）
#     - 或至少标记出来，单独评估

第二、三步：表示与划分

from rdkit.Chem import rdFingerprintGenerator, Descriptors
from rdkit.Chem.Scaffolds import MurckoScaffold
from collections import defaultdict

# ---- 表示：ECFP + 描述符（标准基线）----
gen = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=2048)
DESC = ["MolWt", "MolLogP", "TPSA", "NumHDonors", "NumHAcceptors",
        "NumRotatableBonds", "RingCount", "FractionCSP3",
        "NumAromaticRings", "HeavyAtomCount"]

def featurize(smiles_list):
    X = []
    for s in smiles_list:
        m = Chem.MolFromSmiles(s)
        if m is None:
            X.append(np.zeros(2048 + len(DESC)))
            continue
        fp = np.array(gen.GetCountFingerprint(m).ToList())
        d = np.array([getattr(Descriptors, n)(m) for n in DESC])
        X.append(np.concatenate([fp, np.nan_to_num(d)]))
    return np.array(X)

# ---- 划分：【骨架划分，不要随机划分】（见 051）----
def scaffold_split(smiles_list, frac_train=0.8):
    groups = defaultdict(list)
    for i, s in enumerate(smiles_list):
        m = Chem.MolFromSmiles(s)
        if m is None:
            continue
        scaf = MurckoScaffold.MurckoScaffoldSmiles(mol=m) or f"acyclic_{i}"
        groups[scaf].append(i)
    order = sorted(groups.values(), key=len, reverse=True)
    n_train = frac_train * len(smiles_list)
    train, test = [], []
    for g in order:
        (train if len(train) + len(g) <= n_train else test).extend(g)
    return train, test

# 【为什么必须用骨架划分】：
#   随机划分下，测试集分子的近似物在训练集中
#   → 模型只需「查表 + 插值」
#   → 【R² 可能虚高 0.2~0.3】（见 051）

第四、五步：训练与评测

import lightgbm as lgb
from scipy import stats

def train_and_evaluate(X_train, y_train, X_test, y_test, n_seeds=5):
    """多种子训练，返回集成预测与不确定性"""
    models = []
    for seed in range(n_seeds):
        m = lgb.LGBMRegressor(
            n_estimators=1000, learning_rate=0.05, num_leaves=31,
            subsample=0.8, colsample_bytree=0.8,
            random_state=seed, verbose=-1)
        m.fit(X_train, y_train)
        models.append(m)

    preds = np.array([m.predict(X_test) for m in models])
    mean_pred, std_pred = preds.mean(axis=0), preds.std(axis=0)

    rmse = np.sqrt(((mean_pred - y_test) ** 2).mean())
    mae = np.abs(mean_pred - y_test).mean()
    r2 = 1 - ((mean_pred - y_test) ** 2).sum() / \
             ((y_test - y_test.mean()) ** 2).sum()
    spearman = stats.spearmanr(y_test, mean_pred).correlation

    print(f"RMSE = {rmse:.3f}")
    print(f"MAE  = {mae:.3f}")
    print(f"R²   = {r2:.3f}")
    print(f"Spearman = {spearman:.3f}")
    return models, mean_pred, std_pred

# 【指标的选择要对应决策】（见 155）：
#   排序候选分子 → Spearman、Kendall
#   预测具体数值 → RMSE、MAE
#   虚拟筛选 → 富集因子
#
# 【一个重要的参照】：
#   实验本身的重复性是多少？
#   → 如果同一分子的重复测定标准差是 0.3 个 log 单位，
#     那么 RMSE = 0.3 就已经是【极限】
#   → 【不要追求超过实验噪声的精度】

第六步：适用域

# 【这一步最常被省略，但对实际使用最重要】

from rdkit import DataStructs

def applicability_domain(query_smiles, train_smiles):
    """用最大 Tanimoto 相似度判断适用域"""
    train_fps = [gen.GetFingerprint(Chem.MolFromSmiles(s))
                 for s in train_smiles if Chem.MolFromSmiles(s)]
    results = []
    for s in query_smiles:
        m = Chem.MolFromSmiles(s)
        if m is None:
            results.append({"max_sim": 0.0, "verdict": "无效"})
            continue
        sims = DataStructs.BulkTanimotoSimilarity(
            gen.GetFingerprint(m), train_fps)
        max_sim = max(sims)
        results.append({
            "max_sim": max_sim,
            "verdict": ("域内" if max_sim > 0.6 else
                        "边缘" if max_sim > 0.4 else "【域外】"),
        })
    return results

# 【模型输出应该包含】：
#   预测值 + 不确定性 + 适用域判断 + 最近邻对比
#   → 【只给一个数字的模型，在实际决策中价值有限】（见 155）

QSAR 的适用边界

需要足够的数据：少于 50 个数据点基本无法建模；100～500 可以做粗略的排序；
只在训练的化学空间内可靠：新骨架的预测不可信（见 051《Scaffold Split》、167《模型外推性》）；
不能外推到训练范围之外：训练集最高 pIC50 是 8，模型基本不会预测出 9.5——先导优化后期尤其要注意；
活性悬崖是固有难点：结构相似但活性差异巨大的分子对，模型必然预测错；
精度受实验噪声限制：不可能比实验的重复性更准；
不解释机理：QSAR 给出统计关联，不告诉你为什么——机理需要结构信息（见 096《Binding Pose》）。

关键要点

数据准备与划分的影响远大于调模型——新手常把时间花错地方；
必须转成对数尺度；删失数据（「> 10 μM」）不能当作精确值；
用骨架划分，随机划分会让 R² 虚高 0.2～0.3；
不要追求超过实验噪声的精度；输出必须配套不确定性与适用域判断。

延伸资源

QSPR：045《QSPR 入门》；分子标准化：046《分子标准化》；Scaffold Split：051《Scaffold Split》；
活性预测模型：155《AI 活性预测模型》；不确定性：166《不确定性估计》；ChEMBL：226《ChEMBL》。