4. 分子描述符计算实践：RDKit与PyBioMed实战

分子描述符，说白了就是把化学结构变成数字。我刚开始做QSAR时，总觉得这一步就是跑个脚本的事。直到有一次，我算出来的描述符跟文献对不上，排查了半天——原来是版本问题，RDKit的某个描述符在旧版本里算的是另一个值。

从那天起，我对描述符计算就多了几分敬畏。今天咱们就亲手把这件事做扎实。

4.1 先用RDKit搞定2D描述符

RDKit是目前最主流的化学信息学工具包。我习惯拿它先算2D描述符，因为这些描述符物理意义明确，计算又快。安装方法很简单，直接用pip：

pip install rdkit-pypi

装好之后，来一段最基础的代码。我会先建一个分子对象，然后逐个调描述符函数：

from rdkit import Chem
from rdkit.Chem import Descriptors

mol = Chem.MolFromSmiles('CC(=O)OC1=CC=CC=C1C(=O)O')

# 常见的几个描述符
mw = Descriptors.MolWt(mol)
logp = Descriptors.MolLogP(mol)
hba = Descriptors.NumHAcceptors(mol)
hbd = Descriptors.NumHDonors(mol)
rot_bonds = Descriptors.NumRotatableBonds(mol)
tpsa = Descriptors.TPSA(mol)

print(f"分子量: {mw:.2f}")
print(f"LogP: {logp:.2f}")
print(f"氢键受体: {hba}")
print(f"氢键供体: {hbd}")
print(f"可旋转键: {rot_bonds}")
print(f"拓扑极性表面积: {tpsa:.2f}")

如果你想一次性算出所有可用的2D描述符，可以这样写：

from rdkit.ML.Descriptors import MoleculeDescriptors

# 获取所有描述符的名字
descriptor_names = [desc[0] for desc in Descriptors._descList]
calculator = MoleculeDescriptors.MolecularDescriptorCalculator(descriptor_names)
descriptor_values = calculator.CalcDescriptors(mol)

# 做成字典方便查看
desc_dict = dict(zip(descriptor_names, descriptor_values))
print(f"共计算了 {len(desc_dict)} 个描述符")
print("前5个：", list(desc_dict.items())[:5])

4.2 PyBioMed：多种描述符一把抓

RDKit只能算2D描述符，但实际工作中，我们经常需要其他类型的描述符。PyBioMed这个包整合了多个工具，能算指纹、分子性质、甚至蛋白质描述符。它底层依赖RDKit和其它库，安装命令：

pip install PyBioMed

我还是用阿司匹林举例，对比一下两种工具的计算结果：

from PyBioMed.PyMolecule import PySmiles

# 初始化
calculator = PySmiles.PySmiles()
smiles = 'CC(=O)OC1=CC=CC=C1C(=O)O'

# 计算多种描述符
results = calculator.GetAllDescriptor(smiles)

# 查看有哪些类别
for key in results.keys():
    print(f"- {key}: {len(results[key])} 个描述符")

# 比如看看分子性质部分
molecular_properties = results['MolecularProperty']
print("\n分子性质:",
      {k: v for k, v in list(molecular_properties.items())[:5]})

输出结果一般是这样的：

- MolecularProperty: 27 个描述符
- Constitutional: 22 个描述符
- Topological: 49 个描述符
- Geometrical: 5 个描述符
- Charge: 66 个描述符
- CPSA: 39 个描述符
- MOEtype: 60 个描述符
- Fingerprint: 208 个描述符

嗯，看到这个数字你就明白了——PyBioMed一口气给你算300多个描述符。说实话，全部丢进模型肯定不行。我一般只挑两到三个类别使用，比如分子性质加指纹。

RDKit vs PyBioMed 核心差异

对比维度	RDKit 2D描述符	PyBioMed 全描述符
计算效率	极高，每秒可处理上千分子	中等，指纹和3D描述符较慢
描述符类别	仅2D理论描述符	含指纹、拓扑、几何、电荷等8大类
适用场景	快速筛选、简单模型	复杂的多样本QSAR建模
依赖深度	独立计算，无需外部工具	依赖RDKit + OpenBabel + CDK

4.3 描述符标准化——容易被忽略的关键步骤

描述符算出来了，但它们的尺度差异很大。分子量可能几百，LogP却只有个位数。你想想看，如果直接把原始值扔给机器学习模型，大数值的描述符会主导模型训练。我刚开始学QSAR时就没注意这点，模型跑完后一看系数，全是分子量的贡献，别的描述符基本没起作用。

标准化的方法很多，我推荐两种最实用的：

4.3.1 Z-score标准化（零均值归一化）

这个方法把每个描述符转为均值为0、标准差为1的分布。适用于数据本身近似正态分布的情况。R用scikit-learn实现很简单：

from sklearn.preprocessing import StandardScaler
import pandas as pd

# 假设我们已经算出了描述符矩阵 X
# X 是形状为 (n_samples, n_features) 的 numpy 数组
df = pd.DataFrame(X, columns=descriptor_names)

# 初始化标准化器
scaler = StandardScaler()
X_scaled = scaler.fit_transform(df)

# 转回DataFrame
df_scaled = pd.DataFrame(X_scaled, columns=descriptor_names)

print("标准化前均值:", df.mean().round(2).tolist()[:5])
print("标准化后均值:", df_scaled.mean().round(2).tolist()[:5])
print("标准化后方差:", df_scaled.var().round(2).tolist()[:5])

4.3.2 Min-Max标准化（归一化到[0,1]）

当你的描述符取值有明确上下界时，用Min-Max更直观。比如一些指纹描述符是0到1之间的概率值，用Z-score反而会丢失信息。

from sklearn.preprocessing import MinMaxScaler

minmax_scaler = MinMaxScaler(feature_range=(0, 1))
X_scaled_mm = minmax_scaler.fit_transform(df)

df_scaled_mm = pd.DataFrame(X_scaled_mm, columns=descriptor_names)
print("标准化后范围: [", df_scaled_mm.min().min(), ",", df_scaled_mm.max().max(), "]")

# 正确做法
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 只用训练集拟合
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)

# 用训练集学到的均值和标准差，转换测试集
X_test_scaled = scaler.transform(X_test)

print("训练集均值:", X_train_scaled.mean(axis=0).round(2)[:5])
print("测试集均值（应接近0）:", X_test_scaled.mean(axis=0).round(2)[:5])

4.4 实战要点速查

我把计算和标准化的要点列成清单，方便你以后参考：

• 计算顺序：先确认分子结构无误（检查是否含有特殊原子），再批量计算描述符。
• 检查异常值：计算完第一件事就是检查是否有NaN或Inf值，尤其是用PyBioMed计算电荷描述符时。
• 标准化时机：先划分训练/测试集，再标准化——这个顺序千万别搞反。
• 保存scaler：Fit完StandardScaler后，用pickle存起来。之后在新数据预测时直接load出来transform，保证一致性。

import pickle

# 训练时保存
with open('descriptor_scaler.pkl', 'wb') as f:
    pickle.dump(scaler, f)

# 预测时加载
with open('descriptor_scaler.pkl', 'rb') as f:
    scaler_loaded = pickle.load(f)

X_new_scaled = scaler_loaded.transform(X_new)