第 3 章：分子描述符计算基础——SMILES 与常用描述符库

做 QSAR 的人都知道一句话：垃圾进，垃圾出。描述符算不准，后面模型再花哨也是白搭。这一章咱们就聊聊分子描述符怎么来、怎么算、怎么避开那些坑。我个人习惯把描述符计算分成两步：先把分子结构洗干净，再拿靠谱的工具去算。下面一个个说。

3.1 SMILES：分子的“身份证号”

SMILES 全称是 Simplified Molecular Input Line Entry System。说白了，就是用一串 ASCII 字符把分子结构写出来。你想想看，一个苯环写成一串 "c1ccccc1"，比画个六边形方便太多了。

举个简单的例子：

乙醇：CCO
苯酚：c1ccccc1O
阿司匹林：CC(=O)Oc1ccccc1C(=O)O

SMILES 有几个核心规则：

• 原子用元素符号表示，比如 C、O、N、Cl。芳香性原子用小写，c、o、n。
• 单键、双键、三键分别用 -、=、# 表示，但单键通常省略。
• 分支用圆括号，比如 CC(C)C 表示异丁烷。
• 环结构用数字标记，连接成环的两个原子用相同数字。

3.2 分子描述符：把结构变成数字

描述符是什么？就是把分子的结构信息编码成数值特征。分子量、LogP、氢键供体数……这些都是描述符。我见过有人一口气算 5000 多个描述符，结果模型过拟合得一塌糊涂。所以你想想看，选描述符其实和选特征一样，不是越多越好，是对路才好。

描述符一般分这么几类：

嗯，这里要注意：3D 描述符依赖构象。你算出来的体积是能量最低构象的，还是某个随机构象的？这个差异有时候能大到让你怀疑人生。我在项目里吃过这个亏，后面会细说。

3.3 RDKit：最顺手的描述符工具箱

RDKit 是开源界做化学信息学最成熟的库。没有之一。不管是 SMILES 解析、分子可视化、还是描述符计算，它都是我的首选。

3.3.1 安装 RDKit

我最推荐用 conda 装，省心：

conda install -c conda-forge rdkit

用 pip 也可以，但注意版本依赖：

pip install rdkit-pypi

3.3.2 核心操作：SMILES 转分子对象

from rdkit import Chem

smiles = "c1ccccc1O"
mol = Chem.MolFromSmiles(smiles)

if mol is not None:
    print("分子解析成功")
    print("原子数:", mol.GetNumAtoms())
    print("化学式:", Chem.rdMolDescriptors.CalcMolFormula(mol))
else:
    print("SMILES 无效，请检查")

输出示例：

分子解析成功
原子数: 7
化学式: C6H6O

3.3.3 计算常见的 2D 描述符

RDKit 的 DescriptorCalculator 可以一口气算几百个描述符。我一般先算几个关键的看看：

from rdkit.Chem import Descriptors

mol = Chem.MolFromSmiles("CC(=O)Oc1ccccc1C(=O)O")

mw = Descriptors.MolWt(mol)            # 分子量
logp = Descriptors.MolLogP(mol)        # LogP
hbd = Descriptors.NumHDonors(mol)      # 氢键供体数
hba = Descriptors.NumHAcceptors(mol)   # 氢键受体数
tpsa = Descriptors.TPSA(mol)           # 拓扑极性表面积

print(f"MW: {mw:.2f}, LogP: {logp:.2f}, HBD: {hbd}, HBA: {hba}, TPSA: {tpsa:.2f}")

MW: 180.16, LogP: 1.19, HBD: 1, HBA: 4, TPSA: 63.60

3.3.4 批量计算描述符

实际项目中你不会只算一个分子。我习惯准备一个 CSV，第一列是 SMILES，然后批量算：

import pandas as pd
from rdkit import Chem
from rdkit.Chem import Descriptors

df = pd.read_csv("molecules.csv")
smiles_list = df["smiles"].tolist()

desc_list = []
for s in smiles_list:
    mol = Chem.MolFromSmiles(s)
    if mol is None:
        desc_list.append(None)
        continue
    desc = {
        "MW": Descriptors.MolWt(mol),
        "LogP": Descriptors.MolLogP(mol),
        "TPSA": Descriptors.TPSA(mol),
        "NumHDonors": Descriptors.NumHDonors(mol),
        "NumHAcceptors": Descriptors.NumHAcceptors(mol),
        "NumRotatableBonds": Descriptors.NumRotatableBonds(mol),
    }
    desc_list.append(desc)

df_desc = pd.DataFrame(desc_list)
df_desc.to_csv("descriptors_rdkit.csv", index=False)

3.4 PaDEL-Descriptor：补充 1D 和 2D 描述符的利器

RDKit 算描述符很全面，但有些专业描述符（比如 PubChem 指纹、E-state 指数）它不直接支持。这时候 PaDEL-Descriptor 就派上用场了。它是一个 Java 写的工具，算描述符非常稳。

3.4.1 安装 PaDEL-Descriptor

先去 GitHub 下载 jar 包（搜索 PaDEL-Descriptor）：

wget https://github.com/datapublicate/PaDEL-Descriptor/raw/master/PaDEL-Descriptor.jar

需要 Java 8 以上环境：

java -version  # 确认 Java 已安装

3.4.2 命令行使用

java -jar PaDEL-Descriptor.jar \
  -dir input_sdf/ \
  -file output.csv \
  -descriptortypes all \
  -fingerprints all

参数说明：

• -dir：输入目录，放 SDF 或 MOL 文件
• -file：输出 CSV 文件
• -descriptortypes all：算所有 1D/2D 描述符
• -fingerprints all：算所有指纹类型

也可以用 SMILES 文件直接当输入：

java -jar PaDEL-Descriptor.jar \
  -file input_smiles.smi \
  -output output.csv \
  -descriptortypes all

3.4.3 在 Python 中调用 PaDEL

我通常用 subprocess 调用 jar：

import subprocess
import os

jar_path = "PaDEL-Descriptor.jar"
input_smi = "data.smi"
output_csv = "padel_descriptors.csv"

cmd = [
    "java", "-jar", jar_path,
    "-file", input_smi,
    "-output", output_csv,
    "-descriptortypes", "all",
    "-fingerprints", "all"
]

subprocess.run(cmd, check=True)
print(f"PaDEL 计算完成，结果保存在 {output_csv}")

3.5 本章核心流程总览

来，我们把从 SMILES 到描述符的完整流程捋一遍。下面这张图是我自己项目里常用的工作流：

这张图里说白了就三件事：输入校验、分工具计算、合并输出。我个人建议新手先用 RDKit 算基础描述符，搞不定的再上 PaDEL 补。别一开始就所有工具全上，容易把自己搞晕。

3.6 避坑指南

最后分享几个我在项目中真实遇到的坑：

• SMILES 中的芳香性标记：苯环写成 "c1ccccc1" 没问题，但有些老数据用的是 "C1=CC=CC=C1"。RDKit 能处理，但算出来的指纹可能会有差异。我建议统一用芳香性小写写法。
• 3D 描述符的构象依赖：同样的分子，用 RDKit 的 ETKDG 生成构象和用 OpenBabel 生成，算出来的体积可能差 5%~10%。如果你要用 3D 描述符，务必在论文里写明构象生成方法。
• PaDEL 的内存问题：分子数量超过 10 万时，PaDEL 默认内存可能不够。加参数 -Xmx8g 分配更多内存：java -Xmx8g -jar PaDEL-Descriptor.jar ...
• 描述符的 NA 值处理：有些描述符对特定结构算不出来（比如某些金属配合物）。别直接填 0，建议用均值填充或者删除该列。我试过直接填 0，结果模型学出了一个错误的偏置。