第4章 分子描述符计算：2D描述符实战

讲ADMET建模，绕不开一个基础问题——怎么把分子结构转成数字。

说白了，计算机看不懂结构式，它只认识数字。分子描述符就是那座桥。今天咱们重点聊2D描述符，特别是那三个最核心的：脂溶性、氢键供/受体、分子量。实战中用RDKit和PaDEL两个工具轮着来。

我个人习惯，做小分子ADMET建模，先用RDKit快速算一轮基础描述符，碰到模型效果不好时，再上PaDEL算个大的。后面会细说为什么这么干。

4.1 你需要了解的描述符基础

2D描述符只依赖分子的拓扑结构，不涉及3D构象。这意味着计算快，可重复性好。你想想看，同一个分子在不同力场下3D构象可能不一样，但2D描述符永远是稳定的。这在工业级建模里特别重要。

先列一下咱们这章要重点算的描述符：

描述符 缩写 意义 在ADMET建模中的作用
脂溶性 LogP 正辛醇/水分配系数的对数值 通透性、代谢稳定性、毒性预测
分子量 MolWt 分子质量（Da） 类药性过滤、Caco-2通透性
氢键供体数 HBD 可提供氢原子的个数 血脑屏障渗透、口服生物利用度
氢键受体数 HBA 可接受氢原子的个数 与HBD搭配，影响结合自由能
可旋转键数 RotBonds 非环且非末端可旋转键数 构象灵活度、结合位点适应性
拓扑极性表面积 TPSA 极性原子贡献的表面积 通透性、P-gp底物预测

嗯，这里要注意。LogP有计算值和实验值之分。做建模时，我们一般用计算LogP（cLogP或AlogP），因为它对所有分子都有值。实验中测定的LogP虽然准，但很多化合物根本测不了——比如溶解度太低的。

4.2 用RDKit算2D描述符

RDKit是开源界的瑞士军刀。我最早用RDKit做描述符计算时还踩过坑——它默认的LogP模型是哪一种？后来翻文档才知道，RDKit的cLogP是基于atom‑typing的Wildman‑Crippen方法。这个方法对常见骨架挺好的，但对含氟、含硅的奇葩分子就不太准了。

代码很简单：

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

mol = Chem.MolFromSmiles('CC(=O)Nc1ccc(O)cc1')
if mol:
    logp = Descriptors.MolLogP(mol)
    mw = Descriptors.MolWt(mol)
    hbd = Lipinski.NumHDonors(mol)
    hba = Lipinski.NumHAcceptors(mol)
    rotb = Lipinski.NumRotatableBonds(mol)
    tpsa = Descriptors.TPSA(mol)
    print(f'LogP: {logp:.2f}')
    print(f'MolWt: {mw:.2f}')
    print(f'HBD: {hbd} | HBA: {hba}')
    print(f'RotBonds: {rotb}')
    print(f'TPSA: {tpsa:.2f}')

输出结果：

LogP: 1.10
MolWt: 151.16
HBD: 2
HBA: 2
RotBonds: 1
TPSA: 49.33

我曾经在一次建模项目中，发现模型对含嘧啶环的分子预测偏差很大。查来查去，问题出在RDKit对含氮杂环的LogP计算偏保守。那一次我把数据分成含杂环和不含杂环两组分别建模，效果才上来。这件事让我养成了一个习惯——先算描述符，跑个分布直方图看看，有没有奇怪的离群点。

4.3 PaDEL描述符计算

PaDEL是一个专门的描述符计算工具箱。它支持1D、2D和指纹，总数超过1800个。RDKit算不了的某些分子拓扑特征——比如EState数、AUTOCORR、WHIM描述符——PaDEL全都有。

PaDEL有两种用法，GUI和命令行。做批量计算我推荐命令行。

下载PaDEL-Descriptor.jar后，在终端里跑：

java -jar PaDEL-Descriptor.jar \
    -dir ./input_sdf \
    -file ./output_descriptors.csv \
    -descriptortypes ./PA_2Donly.txt \
    -retainorder \
    -maxruntime 60000

参数说明：

  • -dir：输入的SDF或SMILES文件夹
  • -file：输出CSV路径
  • -descriptortypes：描述符类型文件，控制只算哪些
  • -retainorder：保持输入顺序输出
  • -maxruntime：单个分子超时限制（毫秒）

PA_2Donly.txt文件内容示例：

# 2D描述符类型列表
Aromatic atoms count
Autocorrelation (2D)
Barysz matrix
CPSA (2D)
EState
Extended Topochemical Atom (ETA) indices
FMF descriptor
Fragment Complexity
Ghose-Crippen LogP
HBondAcceptorCount
HBondDonorCount
Molecular weight
Number of Rotatable Bonds
Topological polar surface area
VABC index
WeightedPath
XLogP

输出CSV长这样：

Name,nAcid,nBase,nAtom,ALogP,ALogp2,AMR,Apol,Aroma...
CC(=O)Nc1ccc(O)cc1,0,1,11,1.10,1.21,38.45,15.22,6,...
CC(C)Cc1ccc(C(O)=O)cc1,1,0,13,2.85,8.12,47.33,19.64,7,...

这里有个坑——PaDEL默认把所有描述符全算一遍，近2000列。对机器学习来说，维度灾难不是闹着玩的。我曾经有一批1000个分子，用完PaDEL全算完后，特征矩阵5万列。后来模型只跑了5分钟就过拟合了。所以我强烈建议，先通过描述符类型文件过滤，只选你需要的那几类。

4.4 两个工具怎么选？

我一般这样判断：

  • 快速原型或小规模数据（< 5000分子）：RDKit就够了，代码干净、依赖少、结果好解释。
  • 大规模筛选或需要特殊描述符：PaDEL有绝对优势，比如EState指数在预测hERG毒性时就比RDKit的2D描述符效果好。
  • 生产流水线：我习惯两个都算。RDKit算基础描述符用于基线模型，PaDEL算扩展描述符用于集成模型。
对比维度 RDKit PaDEL
开源/许可证 BSD GPLv2
描述符总数 约200个 约1800个
计算速度 快 中等
接口易用性 Python原生 命令行/Java
指纹支持 有（ECFP等） 有（CDK指纹等）
适合场景 日常建模、特征选择 大规模、细粒度特征工程

4.5 实操锦囊：写一个可复用的描述符计算函数

项目做多了你会发现，描述符计算的逻辑在不同任务间高度相似。我习惯封装成函数，放在一个utils.py模块里反复用。

def compute_rdkit_2d(smiles_list, names=None):
    """
    批量计算RDKit 2D描述符
    """
    from rdkit import Chem, RDLogger
    from rdkit.Chem import Descriptors, Lipinski
    import pandas as pd

    RDLogger.logger().setLevel(RDLogger.ERROR)

    results = []
    for i, smi in enumerate(smiles_list):
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            results.append({f'RDKit_{i}': None for i in range(6)})
            continue
        row = [
            Descriptors.MolLogP(mol),
            Descriptors.MolWt(mol),
            Lipinski.NumHDonors(mol),
            Lipinski.NumHAcceptors(mol),
            Lipinski.NumRotatableBonds(mol),
            Descriptors.TPSA(mol)
        ]
        results.append(row)

    cols = ['LogP', 'MolWt', 'HBD', 'HBA', 'RotBonds', 'TPSA']
    df = pd.DataFrame(results, columns=cols)
    if names is not None:
        df.insert(0, 'Name', names)
    return df

调用：

smiles = ['CC(=O)Nc1ccc(O)cc1', 'CC(C)Cc1ccc(C(O)=O)cc1']
names = ['paracetamol', 'ibuprofen']
df = compute_rdkit_2d(smiles, names)
print(df)

输出：

       Name  LogP   MolWt  HBD  HBA  RotBonds   TPSA
0  paracetamol  1.10  151.16    2    2         1  49.33
1   ibuprofen  2.85  206.28    1    2         3  37.30

你看，两行代码就能把分子转成结构化特征。整章下来就一个核心思路——描述符是建模的起点，起点不对后面全白费。RDKit和PaDEL各有脾性，摸熟了就知道什么时候用谁。

最后给个小提醒：如果你拿到的SMILES有手性信息和同位素标记（比如[13C]），RDKit会无警告地忽略它们，直接按默认同位素算。PaDEL也一样。所以清洗SMILES时，建议先去同位素、去手性，否则同一化合物不同来源的SMILES会算出来不一样的描述符。我就是这么被坑过一次。

← 上一章 📖 返回目录 下一章 →