📚 课程目录第3章分子描述符计算

分子描述符计算：RDKit、Mordred 等工具实战归纳

做 ADMET 预测，说白了，输入特征就是分子描述符。

你模型再花哨，特征没算对，全白搭。我记得刚开始做这个方向的时候，拿 ChemAxon 算了一大堆二维描述符，结果跟 RDKit 对不上，找了三天 bug，最后发现是芳香性判断标准不一样。这种坑，掉过一次就记住了。

为什么要有专门的工具来做这件事？

你想想看，SMILES 字符串只是一串文本，计算机不认识「亲脂性」「氢键供体」这些概念。分子描述符计算工具的本质，就是把化学结构翻译成机器能理解的数值矩阵。

我个人的习惯是：

基础描述符用 RDKit，开源、稳定、文档全。
想用更前沿、更丰富的特征时，再上 Mordred。
如果团队有商业预算，可以考虑 Dragon，但我觉得开源路线完全够用。

一个我踩过的坑： RDKit 的版本差异会导致部分描述符值出现微小变化。我在项目中曾经遇到过 2019 版和 2021 版算出的 LogP 差了 0.02，虽然不大，但在 SVR 调优阶段，这种噪声足以让网格搜索结果不稳定。建议锁定一个版本。

RDKit：安身立命的标配工具

RDKit 是我最常用的工具，没有之一。它的描述符计算模块非常成熟，几乎涵盖了传统 QSAR 需要的所有基本特征。

基础用法

先把环境搭好。我用的是 conda 管理，省心：

conda create -n admet_py310 python=3.10
conda activate admet_py310
conda install -c conda-forge rdkit

然后就能在 Python 里直接调用。我一般会封装一个简单函数：

from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski

def calc_rdkit_descriptors(smiles):
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    return {
        'MolWt': Descriptors.MolWt(mol),
        'LogP': Descriptors.MolLogP(mol),
        'NumHDonors': Lipinski.NumHDonors(mol),
        'NumHAcceptors': Lipinski.NumHAcceptors(mol),
        'TPSA': Descriptors.TPSA(mol)
    }

这样几行代码，就能拿到最基本的五个描述符。但我建议大家不要止步于此，RDKit 内置了超过 200 种描述符。

批量计算技巧

如果手头有几千个分子，一个个循环会慢到怀疑人生。我的做法是利用 Pandas 的 apply 搭配并行：

import pandas as pd
from multiprocessing import Pool

def calc_batch(smiles_list):
    with Pool(processes=8) as pool:
        results = pool.map(calc_rdkit_descriptors, smiles_list)
    return pd.DataFrame(results)

注意： 多进程下 RDKit 的 MolFromSmiles 是线程安全的，但如果你在进程间共享了 molecule 对象，可能会出问题。我最开始就犯过这种错，记住了：每个进程各自创建 molecule，不要传递。

Mordred：新一代描述符计算库

RDKit 虽好，但有些现代描述符（比如 WHIM、GETAWAY）它不直接提供。这时候 Mordred 就派上用场了。

Mordred 的特点是：

纯 Python 实现，安装简单
支持 1825 种描述符（你没看错，将近两千个）
输出格式跟 Pandas 无缝对接

安装同样简单：

pip install mordred

用法跟 RDKit 类似，但 Mordred 返回的是描述符对象，需要转成数值：

from mordred import Calculator, descriptors
from rdkit import Chem

calc = Calculator(descriptors, ignore_3d=False)
mol = Chem.MolFromSmiles('CCCC')
result = calc(mol)

# 转成字典
desc_dict = {str(d): result[d] for d in result}

    我的经验： Mordred 的 3D 描述符需要分子的三维构象，计算会慢很多。在实际的 SVR 调优中，我通常先只用 2D 描述符跑一轮基线，只有预训练表现不错时，再考虑加入 3D 特征。别一上来就把所有描述符全打开，否则你会卡在计算阶段一整天。

两个工具的选择对比

维度	RDKit	Mordred
描述符数量	约 200+	约 1825
计算速度	快（C++ 后端）	较慢（纯 Python）
文档质量	优秀，社区活跃	良好，但示例较少
3D 描述符	需额外模块	内置支持
适用场景	日常快速验证	精细化特征工程

我现在的习惯是：用 RDKit 快速验证想法，用 Mordred 做最终模型的特征工程。两者互补，不冲突。

你可能遇到的两个大坑

嗯，这里要注意。实际操作中，下面两个问题几乎每个人都遇到过。

分子解析失败：SMILES 写错了，或者包含 RDKit 不支持的原子类型（比如铍、锑）。我曾经在处理一组金属有机框架分子时，花了一下午调试——其实是 SMILES 里漏了一个方括号。
描述符全是 NaN：常见于氢原子数、键级相关的描述符。Mordred 在某些情况下会返回缺失值。解决方案是对每个分子做校验，并设置默认填充值。

我会这样处理：

def safe_calc(smiles):
    try:
        mol = Chem.MolFromSmiles(smiles)
        if mol is None:
            return {}
        # 计算所有描述符
        vals = calc(mol)
        # 将 NaN 替换为 0
        return {k: v if v == v else 0.0 for k, v in vals.items()}
    except Exception:
        return {}

一个小技巧： 如果你的分子很多，建议先用 RDKit 的 SanitizeMol 做一遍合法性校验，而不是直接在 Mordred 里报错后再回头检查。这能省下至少一半的调试时间。

描述符计算的完整流程（单章 SVG）

下面这张图展示了我个人在实践中总结的描述符计算与 SVR 调优的衔接流程。你可以把它当作一个对照清单，每次做实验前走一遍。

总结

我个人的建议很明确：新手先用 RDKit 跑通整个 pipeline，等你把 SVR 调优的逻辑搞清楚了，再回头用 Mordred 扩展特征集。

计算描述符看起来是体力活，但这一步要是出了偏差，后面的调优就是缘木求鱼。花点时间把这一步做扎实，比后面盲目调参要有效得多。