2.4 分子描述符：把可测性质压缩成数字

2.4.1 描述符的定义：把分子写成数

2.1 节把分子存成图，2.2 与 2.3 节让它能写、能查。机器学习还差一步：模型吃的是数，不是图。分子描述符(molecular descriptor)补上这一步——本章讨论的二维描述符是一个确定性映射 f，把分子图 G 变成实数 d = f(G)，或一列实数（向量）。“确定性”的准确含义是：算法与参数固定、输入固定，任何时刻重算数值不变；输入是经 sanitize 的分子图时，原子的编号无论怎么排，数值也不变——这正是 1.2 节讨论过的置换不变性(permutation invariance)，在那里它是对模型结构的要求，在这里它先落在特征的定义上。还要预留一个口子：描述符不止二维一族——RDKit 的 Descriptors3D（主转动惯量 PMI、回转半径、偏离球形度等）要先生成三维构象，d = f(G, X) 依赖构象坐标 X，只保证对整体平移与旋转不变；同一分子图的两个构象可以给出不同的数。“描述符是分子图的函数”这个等式，只对二维成立。

并非任何映射都配称好描述符。评价标准可以概括为三条，但先说明它们是挑选代表的偏好，不是定义的门槛。其一，可复算：算法与参数固定，不同实验室、不同软件版本重算应当一致，至少同一实现内严格一致。其二，有化学含义：数值的升降对应可指认的结构变化——分子量涨是因为多了原子，TPSA 涨是因为添了极性基团；含义间接的数学指标（信息论一类）解释与排错都更难，仍是合法描述符，本节挑选代表时优先取含义清晰者。其三，对结构微变敏感而稳定：这是对连续描述符的期望——骨架上换一个取代基，数值应当移动但不抖动；离散计数类（供体数、环数）本来就按整数跳变，不在此列。三条同时做到顶格的描述符并不存在，取舍贯穿整节。

2.4.2 组成类：直接可解释的计数

第一类描述符只做计数与求和：分子量、重原子数、N/O 占比、氢键供体(hydrogen-bond donor)与氢键受体(hydrogen-bond acceptor)计数。化学含义写在定义里。分子量是全原子按同位素自然丰度加权的质量，回答“多大”；重原子数给出规模的粗粒度；N/O 占比 (N+O)/重原子数是极性的粗画像（RDKit 的 NumHeteroatoms 口径更宽，数所有非碳非氢原子）；供体与受体计数回答“能形成几根氢键”——但口径要分清：现代供体计数（RDKit CalcNumHBD）数的是供体原子，一个 NH₂ 算一个位点，而不是两根 N–H；历史 Lipinski 口径（CalcNumLipinskiHBD）才按 N–H/O–H 键数。受体的化学定义是“孤对可供电子的原子”，也不限于氮氧——RDKit 的受体 SMARTS 覆盖 O、S、N 与芳香 n/o/s。

代入咖啡因 C8H10N4O2：分子量 194.19，重原子 14，N/O 占比 6/14 ≈ 0.43，供体 0——三个 N–CH₃ 封掉了全部 N–H。受体的口径差在咖啡因上最典型：历史 Lipinski 计数（CalcNumLipinskiHBA / Lipinski.NOCount）直接数氮氧，得 6；现代计数（CalcNumHBA，Descriptors.NumHAcceptors 与 Lipinski.NumHAcceptors 都走它——后一个名字极具误导性）按受体环境判定，两个酰胺型 N–CH₃ 不算，得 3（两个羰基氧加一个咪唑型氮）。Lipinski “五规则”用的是 HBA ≤ 10 这类粗边界：MW ≤ 500、LogP ≤ 5、HBD ≤ 5、HBA ≤ 10。计数描述符看似最客观，其实口径一变数字就变，这一伏笔留到 2.4.6。三类代表连同机理与局限，先列成表 2.4-1，随后三小节逐一展开。

下面一段代码把表中的四列变成实际输出，顺带印证定义里的“确定性”——重算千遍，数值不动。

from rdkit import Chem
from rdkit.Chem import Descriptors

caf = Chem.MolFromSmiles("CN1C=NC2=C1C(=O)N(C)C(=O)N2C")
funcs = {"MolWt": Descriptors.MolWt, "TPSA": Descriptors.TPSA,
         "MolLogP": Descriptors.MolLogP, "BalabanJ": Descriptors.BalabanJ}
for name, fn in funcs.items():
    print(f"{name:8s}{fn(caf):9.2f}")   # 194.19 / 61.82 / -1.03 / 3.22

2.4.3 物理化学估计类：加和模型

第二类描述符估计一个可测的物理化学量，机理是加和模型(additive model)：把整体性质拆成部件贡献之和，部件是原子或片段，贡献值在实验数据集上回归拟合得到。两个代表：MolLogP 与 TPSA。

MolLogP：Wildman–Crippen 原子贡献法。LogP（分配系数对数）度量分子在正辛醇与水两相间的分配偏好，是亲脂性的标准度量，也是 Hansch 分析以来定量构效关系里出场最多的变量。该法把每个原子——重原子与氢都算——归入数十类原子类型(atom type)之一：类型由元素、杂化方式与连接环境共同决定，如“芳香碳，带一个氢”与“sp³ 碳，连一个氮”是两类；氢同样分型——连碳的氢、连氮的氢、连氧的氢各有系数，并非一个统一的氢贡献（Wildman & Crippen, 1999）。

图 2.4-2 把咖啡因的账目逐原子摊开：三个 N–CH₃ 的碳各 −0.144，四个环氮各 −0.324，两个羰基碳各 −0.819，两个羰基氧各 +0.179，两个稠合芳香碳各 +0.295，C8 的芳香 CH +0.158；重原子小计 −2.259，加 10 个 H1 型氢的 +1.230，得 −1.03。实验值约 −0.07，偏差接近一个对数单位——这是加和模型的系统误差，成因见 2.4.5。

TPSA：Ertl 片段贡献法。拓扑极性表面积(topological polar surface area, TPSA)的原型要先生成三维构象、再算极性原子（氮、氧及其相连氢）的范德华表面积。Ertl 的办法是绕开构象：以极性原子为中心定义数十类片段，在药物分子三维 PSA 的基础上拟合每类片段的面积贡献，直接在拓扑上求和（Ertl et al., 2000）。RDKit 的默认实现只计氮氧片段；Ertl 原始参数还覆盖硫与磷，须 includeSandP=True 才纳入。经验上，不同研究给出过 120–140 Ų 量级的口服暴露界限：Kelder 等观察到经被动跨细胞途径吸收的口服药物 PSA 多不超过约 120 Ų；Veber 等则把 PSA ≤ 140 Ų（连同旋转键 ≤ 10）与较好的口服生物利用度统计关联。这些是数据集依赖的经验规则，不是普适的膜扩散临界值。代价同样来自拓扑——构象把极性原子折叠进分子内部时，真实暴露面积变小，TPSA 对此一概不知。

可旋转键。可旋转键(rotatable bond)的基础定义有三道过滤：单键、不在环上、两端各接非末端重原子；RDKit 默认的严格口径再排除酯、酰胺等（乙酸甲酯的酯 C–O 在严格口径下不计，宽松口径 NonStrict 才计）。它作分子柔性的 proxy。咖啡因为 0：稠合双环把所有键锁死。柔性影响结合时的熵代价与选择性，药物化学里“刚性化”是常见优化思路，这个计数是它的量化注脚。

2.4.4 拓扑与图论类：图不变量

第三类描述符从连接方式取数。“拓扑”指的是不依赖三维坐标、只在图上计算，而不是“不看元素”——本节要讲的三种指数里，只有最纯粹的拓扑指数只看无标号的邻接关系；χ 的价版本按元素取价电子数，κ 有杂原子修正项，Balaban J 在不饱和体系里用键级加权距离。它们共同的底线是成为图不变量(graph invariant)：图的画法、原子的编号都不影响数值，只有结构本身起作用。分子图的邻接与距离关系经代数变换压成单个数，即拓扑指数(topological index)。

Balaban J 指数。先由距离矩阵求每个原子的距离和(distance sum)s_i（到其余各原子距离之和），再对每条边取 (s_is_j)^−1/2 求和，乘 m/(μ+1) 归一，μ 为环秩(cyclomatic number)，即独立环数（μ = m−n+1 对连通图成立；多组分盐按 μ = m−n+c 累加）。距离怎么算有一条关键细节：饱和烷烃里每条键的长度都是 1；对不饱和与芳香体系，Balaban 原文与 RDKit 的实现都把键级计入——每条键按键级的倒数计长（芳香键 1.5 记 2/3，双键记 1/2），于是同一张图的距离和更小、J 更大。

χ 连接性指数与 κ 形状指数。分子连接性指数(molecular connectivity index)χ 把“分支度”参数化：对指定长度的路径求 (Πδ)^−1/2 之和，δ 为顶点度；价版本 δ^v 用价电子数减氢数替代度，杂原子与重键于是有了区分。

形状指数(shape index)κ 走另一条路：数特定长度的路径。以全单键的烷烃骨架为例，κ₁ = n(n−1)²/P₁²，P₁ 是一阶路径（边）数；无环图恒有 P₁ = n−1，于是 κ₁ 恒等于 n，只测大小。κ₂ = (n−1)(n−2)²/P₂² 才分辨形状：丁烷 P₂ = 2，κ₂ = 3.0；异丁烷 P₂ = 3，κ₂ = 1.33。越伸展数值越大，越缩拢越小（Kier & Hall, 1986；多重键与杂原子另需修正项）。

这一家族的真正卖点，是补上加和类的盲区。表 2.4-2 里茶碱、可可碱、副黄嘌呤互为二甲基位置异构体：组成相同、原子类型相同，只是位置重排。表中这些加和类描述符塌缩成同一组数——加和模型并非原则上分不开位置异构体，片段类型若编码足够大的局部环境也能区分，但常见的原子型加和确实分不开——唯 J 指数把三者分开。组成与加和回答“有什么”，拓扑回答“怎么连”。

再看“微改则微变”。茶碱到咖啡因只多一个 N–CH₃：分子量 +14.03（一个 CH₂），供体从 1 到 0，TPSA 降 10.9 Ų（习题 2.4-1 的机理），MolLogP 几乎不动——甲基的亲脂贡献与 N–H 换 N–CH₃ 的类型变化恰好抵消。一组描述符同时给出四个方向不同的读数，这正是多描述符并用的价值。下面一段代码复算表 2.4-2。

xanthines = [("咖啡因 caffeine 1,3,7-Me", "CN1C=NC2=C1C(=O)N(C)C(=O)N2C"),
             ("茶碱 theophylline 1,3-Me", "CN1C2=C(C(=O)N(C1=O)C)NC=N2"),
             ("可可碱 theobromine 3,7-Me", "CN1C=NC2=C1C(=O)NC(=O)N2C"),
             ("副黄嘌呤 paraxanthine 1,7-Me", "CN1C=NC2=C1C(=O)N(C(=O)N2)C")]
for name, smi in xanthines:            # 同一黄嘌呤骨架，仅甲基位置不同
    m = Chem.MolFromSmiles(smi)
    print(f"{name}: MW={Descriptors.MolWt(m):.2f} LogP={Descriptors.MolLogP(m):.2f}"
          f" TPSA={Descriptors.TPSA(m):.1f} J={Descriptors.BalabanJ(m):.3f}")

输出里三个二甲基异构体（茶碱、可可碱、副黄嘌呤）的 MolWt、LogP、TPSA 一字不差，J 却各不相同——表 2.4-2 的分辨率差异由代码直接印证（咖啡因是三甲基取代，作对照组）。局限同样清楚：图不变量没有化学单位，数值不直接对应任何可测量；判别力强不等于因果关系，J 能分开异构体，却不能说明为什么分开。它们是紧凑的结构摘要，作回归输入尚可，作机理解释要谨慎（Todeschini & Consonni, 2000）。

2.4.5 加和模型的边界与表示的谱系

加和模型的假设只有一条：部件贡献彼此独立。化学里处处是反例。远程相互作用：共轭贯通整个环系时，局部原子类型捕捉不了整体电子结构；分子内氢键把极性基团互相“抱住”，分子外观的亲脂性高于各部件之和。构象依赖：logP 与极性表面积的实验值在溶液中测得，分子采取什么构象、极性面暴露多少，加和模型一概不知。电荷形式：两性离子(zwitterion)（如氨基酸内盐）在水与正辛醇中的行为差别悬殊——实验上随 pH 显著变化的是 logD，给定 pH 下全部分子形式合起来的分配比；logP 专指中性形式本身的分配，原则上不随水相 pH 变。加和模型则完全按输入的结构计算：Crippen 参数表里有各种带正电、负电的原子类型，输入什么质子化形式，就得到那个形式的估计——结果怎么解读，取决于预处理把分子统一成了哪种形式。咖啡因近一个对数单位的偏差是模型的残差：原子类型表没有捕捉这个稠合芳香环系的电子环境，属于近似在此处的系统偏差；对单个分子，残差很难精确拆因。判断何时可信的经验法则：片段彼此独立、无强分子内作用、输入形式与目标量匹配时，加和估计可靠；三者任一不满足，估计只能当量级参考。

信息论视角下，描述符是人工特征：把化学先验手工压缩进低维。压缩必然丢信息——表 2.4-2 的异构体折叠就是明证，芳香族的立体化学、构象群更是整块丢失。把它放进表示方法的谱系里看，边界更清楚：

2.4.6 使用中的陷阱：共线性、口径与 p≫n

共线性。描述符不是相互独立的测量。分子量与重原子数几乎共线；摩尔折射率、MolLogP 与分子尺寸的相关在同类系列里也常很高——但跨化学空间会明显下降，不存在普适的相关阈值。一个“拓扑指数族”内部的相关更是天然。线性回归里这叫多重共线性(multicollinearity)：系数在相关列之间随机分摊，训练集一换，符号就能翻转——“某描述符系数为正”的化学解读随之失去根基。信息没有凭空多出来：一百个高度相关的列，信息量未必超过十个。

口径与分子形式。2.4.2 已经埋了伏笔：同一个“氢键受体数”，Lipinski 计数给 6，现代计数给 3——定义不同，两个数都对，数据表必须记录口径（且 API 名会骗人：Lipinski.NumHAcceptors 走的是现代计数）。更深的坑在分子本身：盐形式不同、互变异构体(tautomer)不同，图就不同，描述符随之改变。TPSA 对 N–H 的位置有多敏感，表 2.4-2 已给出——一个 N–CH₃ 与 N–H 的交换移动 10.9 Ų。麻烦在于标签：实验若测的是另一种形式，特征与标签便错位。2.6 节的清洗工序——盐拆分、主互变异构体选择、价态规范化——需要一个一致且与任务匹配的策略：多数建模任务确实要去盐、选主互变异构体；但预测盐型溶解度或特定互变异构体的性质时，抹掉这些差异反而丢掉标签信息。怎么标准化是任务决策，不是无条件的默认流水线。

p≫n。商业与开源软件能给出的描述符以千计，一个药物项目的活性数据常只有几百条。变量多于样本、设计矩阵又无退化列时，普通最小二乘可以完美内插训练集，纯噪声也照拟合不误；训练 R² 接近 1 只说明自由度耗尽，不说明任何预测力。习题 2.4-3 把口径与自由度两条线索合在一起。

下一节走向谱系的另一端：摩根指纹把圆形结构环境展开成高维稀疏向量，Tanimoto 系数为这类向量配备专门的相似度度量。手工低维与手工高维的分野，在那里看得最清楚。