分子描述符计算——从分子结构到数字特征

分子描述符计算，说白了就是把分子结构变成一串数字。这步做不好，后面的模型全是白搭。我在刚接触QSAR时，总觉得描述符越多越好，结果模型过拟合得一塌糊涂。后来才明白，描述符的质量和选择，比数量重要得多。

你想想看，一个分子用SMILES表示也就是一串字符，但计算机不认识什么苯环、羰基。它只认数字。所以我们要用描述符把化学知识翻译成特征向量。这个翻译过程，决定了模型能学到什么。嗯，这里要注意：描述符计算不是简单的代码搬运，你得理解每个描述符的物理化学含义。

什么是分子描述符？

分子描述符是对分子结构进行数学化表征的结果。它可以是分子量、LogP这种单一数值，也可以是分子指纹这样的高维向量。我习惯把描述符分成三类：

• 2D描述符：从分子拓扑结构计算，不依赖3D构象。比如分子量、氢键供体/受体数、TPSA。计算快，稳定性好。
• 3D描述符：需要分子的三维构象。比如分子体积、溶剂可及表面积、3D分子指纹。计算慢，构象依赖性强。
• 指纹描述符：用位向量表示分子子结构的存在与否。比如MACCS指纹、Morgan指纹（ECFP）。适合做相似性搜索和分类模型。

RDKit：我的主力计算工具

RDKit是目前最成熟的开源化学信息学工具包。我个人习惯用Python调用RDKit，把计算流程写成脚本，方便复用和调试。安装很简单：

conda install -c conda-forge rdkit

或者用pip：

pip install rdkit-pypi

装好之后，加载分子只需要一行代码：

from rdkit import Chem
mol = Chem.MolFromSmiles('CC(=O)Oc1ccccc1C(=O)O')

这里有个坑：MolFromSmiles返回None表示SMILES不合法。我建议每次都做空值检查，尤其是从数据库读数据时。

常用描述符一览

RDKit内置了200多个描述符。我用得最多的是下面这些：

动手算一算——代码实战

下面是我在项目中常用的计算模板。先从单个分子开始：

from rdkit import Chem
from rdkit.Chem import Descriptors, rdMolDescriptors

mol = Chem.MolFromSmiles('CC(=O)Oc1ccccc1C(=O)O')
if mol is None:
    raise ValueError('Invalid SMILES')

# 基础2D描述符
mw = Descriptors.MolWt(mol)
logp = Descriptors.MolLogP(mol)
hba = Descriptors.NumHAcceptors(mol)
hbd = Descriptors.NumHDonors(mol)
tpsa = Descriptors.TPSA(mol)
rotb = Descriptors.NumRotatableBonds(mol)

# Morgan指纹 (ECFP4)
fp = rdMolDescriptors.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048)

print(f'MW: {mw:.2f}, LogP: {logp:.2f}, TPSA: {tpsa:.2f}')

批量计算时，我习惯用pandas组织数据：

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

def calc_descriptors(smiles_list):
    results = []
    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            results.append(None)
            continue
        results.append({
            'MolWt': Descriptors.MolWt(mol),
            'LogP': Descriptors.MolLogP(mol),
            'TPSA': Descriptors.TPSA(mol),
            'HBA': Descriptors.NumHAcceptors(mol),
            'HBD': Descriptors.NumHDonors(mol),
            'RotB': Descriptors.NumRotatableBonds(mol)
        })
    return pd.DataFrame(results)

df = calc_descriptors(['CC(=O)Oc1ccccc1C(=O)O', 'c1ccccc1'])
print(df.head())

记得有一次，我帮同事调试一个模型，发现他的描述符矩阵里混了大量NaN。排查了半天，原来是他的SMILES里有几根不合法键。从那以后，我每次批量计算都会加空值检查和异常捕获。

描述符标准化——别让量纲坑了你

描述符的取值范围差异很大。MolWt在几百，而LogP在正负几。如果不做标准化，模型会把量纲大的特征当成重要特征。我建议在计算描述符后立即做标准化：

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
df_scaled = scaler.fit_transform(df)
df_scaled = pd.DataFrame(df_scaled, columns=df.columns)

这里有个细节：标准化参数（均值和标准差）必须从训练集计算，然后应用到测试集。千万别在全集上做标准化，否则会引入数据泄露。我见过不止一个新手在这里翻车。

描述符选择——少即是多

描述符不是越多越好。冗余和噪声特征会拉低模型泛化能力。我常用的筛选方法有：

• 低方差过滤：去掉几乎所有样本取值相同的描述符。这种特征没有区分度。
• 相关性分析：两两描述符相关性超过0.9，只保留其中一个。比如分子量和摩尔折射率通常高度相关。
• 基于模型的选择：用随机森林或Lasso训练一次，根据特征重要性排序，保留Top N。

我在项目中的做法是：先算50-80个2D描述符，加上2048位Morgan指纹。然后用相关性过滤去掉冗余，再用Lasso做自动选择。这样既能保留化学可解释性，又能控制特征维度。

上面这张图是我自己总结的流程框架。你看，从分子结构到QSAR模型，中间最核心的就是描述符计算和特征选择这两个环节。我每次开始一个新项目，都会先跑一遍这个流程，看看哪些描述符稳定、哪些描述符容易出问题。