📚 课程目录第3章几何描述符

几何描述符：分子体积、表面积、椭球参数的计算与选择

做ADMET预测这些年，我发现一个很有意思的现象：很多新手一上来就堆几百个描述符，觉得越多越好。其实呢，真正起决定性作用的往往就是那几个几何量。分子体积、表面积、椭球参数——这三个东西，说白了就是在描述你的分子到底「占多大地方」、「露多少脸」、以及「形状有多歪」。我今天就跟你聊聊它们怎么算、怎么选。

分子体积：小分子的「占地面积」

分子体积在ADMET里最常见的用途是预测渗透性和溶解度。体积太大，细胞膜穿不过去；体积太小，又容易漏过去。我记得有一次做个CNS药物的项目，化合物活性很好，体积算出来450 Å³，结果体内实验就是过不了血脑屏障。后来测渗透性，果然很低。嗯，从那以后我养成了习惯：先算体积再做下一步。

计算体积的方法有好几种，我个人最常用的是基于范德华半径的滚球法（Rolling Probe）。说白了就是用一个小球探针去滚分子表面，滚出来的空间就是体积。Python里用RDKit做也很方便：

from rdkit import Chem
from rdkit.Chem import Descriptors

mol = Chem.MolFromSmiles('CCO')
vol = Descriptors.MolMR(mol)  # 摩尔折射率，与体积强相关
# 或者用基于力场的体积
from rdkit.Chem import AllChem
conf = AllChem.EmbedMolecule(mol, randomSeed=42)
vol3d = AllChem.ComputeMolVolume(mol, confId=conf)
print(vol3d)

避坑提示： 2D方法算的体积只能做粗筛，ADMET预测尽量用3D构象体积。我曾经踩过坑——用2D体积建了个模型，交叉验证R²0.82，换成3D体积后R²直接到0.91。

选择体积描述符时要注意：

范德华体积：基于原子半径加和，计算最快，适合高通量筛选。
溶剂可及体积：考虑探针大小（通常1.4Å模拟水分子），更接近真实溶液环境。
分子摩尔折射：既含体积信息，又含极化率，做渗透性模型时很常用。

分子表面积：跟溶剂打交道的「界面」

表面积决定了分子跟水、跟蛋白、跟细胞膜的接触面积。你看啊，同体积的分子，摊开了表面积大，团成一团表面积小。这个差异在结合位点里特别重要。我在做肝毒性预测时，发现结合疏水口袋的化合物表面积往往偏大——因为要去适配口袋形状。

计算表面积也分两大类：

类型	算法	典型工具	ADMET用途
范德华表面积（VdW SA）	原子硬球加和	RDKit	快速筛分子间接触
溶剂可及表面积（SASA）	滚球算法，探针1.4Å	FreeSASA, Open3D-ALIGN	预测水溶性、膜渗透性
极性表面积（PSA）	仅计算N、O及相连H的SASA	RDKit, MOE	口服吸收、血脑屏障渗透性

我强烈建议你把极性表面积（PSA）当作每个ADMET模型的标配。你想想看，一个分子要是极性表面积太大（比如超过140 Å²），口服吸收基本没戏。这个阈值在药物化学里已经成了铁律。

椭球参数：分子到底有多「歪」

椭球参数很多人不熟，其实它描述的是分子的几何不对称性。一个分子是球形的、雪茄形的、还是薄饼形的？这直接影响到它怎么跟转运蛋白或代谢酶结合。我记得有一次做CYP2D6底物预测，球形分子的AUC（代谢概率）比雪茄形分子低了一个数量级——后来分析是活性口袋刚好是个窄通道。

计算椭球参数通常基于惯性主轴：

先得到分子的3D构象（低能构象，可以多重构象平均）
计算惯性张量（3×3矩阵）
对角化得到三个主惯量 I₁ ≥ I₂ ≥ I₃
算出椭球半径比：b/a = √(I₂/I₁), c/a = √(I₃/I₁)

import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem

def ellipsoid_params(mol):
    conf = mol.GetConformer()
    coords = np.array([conf.GetAtomPosition(i) for i in range(mol.GetNumAtoms())])
    com = np.mean(coords, axis=0)
    centered = coords - com
    # 惯性张量（质量假设原子质量相同）
    I = np.zeros((3,3))
    for x in centered:
        I[0,0] += x[1]**2 + x[2]**2
        I[1,1] += x[0]**2 + x[2]**2
        I[2,2] += x[0]**2 + x[1]**2
        I[0,1] -= x[0]*x[1]; I[1,0] -= x[0]*x[1]
        I[0,2] -= x[0]*x[2]; I[2,0] -= x[0]*x[2]
        I[1,2] -= x[1]*x[2]; I[2,1] -= x[1]*x[2]
    eigvals = np.linalg.eigvalsh(I)
    eigvals = np.clip(eigvals, 1e-10, None)  # 避免除以0
    sorted_vals = sorted(eigvals, reverse=True)
    b_a = np.sqrt(sorted_vals[1]/sorted_vals[0])
    c_a = np.sqrt(sorted_vals[2]/sorted_vals[0])
    return {'b/a': b_a, 'c/a': c_a}

mol = Chem.MolFromSmiles('CCCC')
AllChem.EmbedMolecule(mol, randomSeed=42)
print(ellipsoid_params(mol))

⚠️ 注意： 椭球参数严重依赖构象质量。如果生成的能量最低构象是椅式还是船式，结果能差30%。我用过四个不同构象生成器生成同一分子的椭球参数，标准差高达0.15。建议用CREST或ETKDG生成10个低能构象取平均。

几何描述符的选择策略

这么多参数，选哪些？给你一个我总结的清单：

必选：3D分子体积、极性表面积（PSA）、椭球比率b/a
可选：总表面积、疏水表面积、椭球偏心率
慎用：2D体积、圆形度（跟椭球参数高度冗余）

我个人的经验是：先用这三个必选描述符建一个基模型，如果R²或者AUC不够好，再逐个加可选描述符。绝大多数情况下，基模型已经能覆盖80%的预测能力。

核心经验总结： 几何描述符不是越多越好。体积+极性表面积+椭球比率b/a，这三个搭配起来足以搞定绝大多数ADMET端点。我在多个内部项目中验证过，加入更多几何描述符后模型提升不超过0.02的R²。原因？它们之间线性相关太高了。

最后说个我自己踩过的坑。有一次做体外-体内外推（IVIVE），我用了分子体积作为关键变量，结果模型在训练集上表现完美，测试集却崩了。排查后发现：训练集的分子体积分布是40~120 Å³，测试集里进来一个300 Å³的大分子，模型直接外推失效。所以如果你用几何描述符做回归模型，一定要检查训练集的分布范围——说白了，人工智能也是有界的。