3、分子描述符特征工程：使用RDKit计算分子指纹与理化性质

做ADMET预测，有一句话我特别认同：数据质量决定模型上限，特征工程决定下限。

你可能会用XGBoost调出一手好参数，但如果喂进去的特征一塌糊涂，那结果基本就是白费功夫。我在几个实际项目里都吃过这个亏——一开始急着调参，后来发现是分子描述符选错了。

说白了，描述符就是我们模型的「眼睛」。它看不清楚，模型就不可能看对。这一章我们重点聊聊怎么用RDKit把分子转化成XGBoost能消化的数字。

3.1 为什么需要分子指纹？

先问个问题：一个分子在计算机里是什么？

通常我们拿到的是SMILES或SDF文件，但XGBoost不认识这些字符串。它只认数值。所以必须把分子结构编码成固定长度的向量。

分子指纹就是干这个的。

它是一种位向量（bit vector），每一位代表子结构是否存在。比如ECFP4指纹，它会遍历分子中每个原子的局部环境，生成一系列哈希值，再映射到固定位长上。

核心原则：指纹的「分辨率」要匹配你的预测任务。对于ADMET，1024位或2048位的ECFP4是主流选择。位长太短容易碰撞，太长又浪费内存。

我记得有个案例：某团队用ECFP6做血脑屏障预测，效果始终不理想。后来换成ECFP4，性能直接提升5个点。为什么？因为ECFP6太细致了，捕捉了太多噪声。你想想看，ADMET的机制本身就很粗糙，太细的指纹反而坏事。

3.2 用RDKit计算分子指纹

RDKit是处理分子的瑞士军刀。安装它很简单：

pip install rdkit-pypi

然后我们写个脚本来计算指纹。下面这段代码我经常在项目里直接用：

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

def smiles_to_ecfp(smiles_list, radius=2, n_bits=2048):
    """
    将SMILES列表转换成ECFP指纹矩阵
    radius=2 对应 ECFP4
    """
    fingerprints = []
    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            # 遇到无效SMILES，用全零向量代替
            fingerprints.append(np.zeros(n_bits))
            continue
        fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius, nBits=n_bits)
        arr = np.array(fp)
        fingerprints.append(arr)
    return np.array(fingerprints)
小技巧：Chem.MolFromSmiles 如果遇到无法解析的SMILES会返回 None，一定要加校验。我之前在项目中因为没加这个，跑了三天模型才发现数据少了10%……

除了ECFP，还有几种指纹值得关注：

指纹类型参数适用场景
ECFP (Morgan)radius=2, nBits=1024-2048ADMET通用，最推荐
FCFPradius=2, 使用特征原子侧重理化性质预测
MACCS166位固定指纹快速对比，丢失信息多
Topological Torsion4原子路径指纹立体化学敏感任务

我个人习惯优先选ECFP4。它平衡了覆盖率和特异性。如果算力允许，我会同时尝试radius=1和radius=3做对比实验。

3.3 理化性质：让模型看懂分子「性格」

指纹描述的是「结构」，而理化性质描述的则是「行为」。比如：

  • logP（脂水分配系数）：影响膜透过性和溶解度
  • 分子量：大的分子通常不易吸收
  • 氢键供体/受体数量：影响与靶点结合
  • 可旋转键数量：影响柔性和代谢
  • 极性表面积 (TPSA)：与肠道吸收强相关
  • 电荷、芳香环数量：影响毒性和代谢

RDKit的Descriptors模块可以一口气算200多个性质。但别全用上——相信我，我试过。结果模型过拟合得一塌糊涂。关键在于筛选。

from rdkit.Chem import Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors

def calc_all_descriptors(smiles_list):
    """
    计算一组理化性质，返回DataFrame
    """
    calc = MoleculeDescriptors.MolecularDescriptorCalculator([
        'MolWt', 'MolLogP', 'NumHDonors', 'NumHAcceptors',
        'NumRotatableBonds', 'TPSA', 'FractionCsp3',
        'NumAromaticRings', 'NumAliphaticRings'
    ])
    
    # 注意：描述符顺序要记牢
    desc_names = calc.GetDescriptorNames()
    results = []
    
    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            results.append([0]*len(desc_names))
            continue
        vals = calc.CalcDescriptors(mol)
        results.append(vals)
    
    import pandas as pd
    df = pd.DataFrame(results, columns=desc_names)
    return df
曾经踩过的坑：MolLogP 计算的是野生动物的logP，不是实验室测定的那个。两者有偏差。如果你有实测logP值，要优先用实测值去校正或替代计算值。我曾经在Caco-2渗透性预测里吃了这个亏，模型精度差了10%。

一般来说，筛选理化性质的原则是：

  1. 相关性优先：先查文献，看哪些性质在ADMET任务中被反复提及
  2. 去冗余：比如NumHAcceptors和NumHDonors有时高度相关，只保留一个
  3. 非负约束：不同尺度的性质要标准化，否则XGBoost的树模型会偏向大尺度的特征

3.4 组合特征：指纹 + 理化性质

真正效果好的是组合拳。拿ECFP4当基础，再加上10-15个精选理化性质，结果往往比只用指纹或只用性质好一大截。

具体的组合方式很简单：

import numpy as np

# 假设我们已经有了指纹矩阵 fp_matrix (n_samples, 2048)
# 和理化性质矩阵 prop_matrix (n_samples, 12)

combined = np.hstack([fp_matrix, prop_matrix])
# 输出 shape: (n_samples, 2060)

这里有个细节：要不要对指纹做标准化？

指纹是0/1的位向量，不需要标准化。但理化性质的尺度差异很大，比如MolWt范围在几十到上千，而NumHDonors在0到10之间。我会用StandardScaler处理理化性质部分：

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
prop_scaled = scaler.fit_transform(prop_matrix)

combined_final = np.hstack([fp_matrix, prop_scaled])
为什么只缩放理化性质？因为指纹的每一位是独立的子结构存在性判断，不服从任何分布假设。强行标准化反而会破坏它的语义。

3.5 特征工程的整体流程

下面这张图总结了从原始SMILES到最终特征矩阵的完整链路，我在项目里一直按这个走：

RDKit 分子描述符特征工程流程 输入：SMILES / SDF Chem.MolFromSmiles 解析 计算分子指纹 ECFP4 / MACCS / FCFP 计算理化性质 logP / TPSA / 分子量 指纹筛选 & 位长调整 性质筛选 & 标准化 最终特征矩阵 → XGBoost

3.6 避坑：特征冗余与维度灾难

我做ADMET项目时，一开始总想「全都要」。结果是特征维度上万，模型训练时间翻倍，精度却不升反降。

你想想看，XGBoost虽然自带特征重要性排序，但它不会帮你节省内存。每多一个无用的特征，树的深度和分裂点数量都可能膨胀。

我的经验法则：
✅ 指纹位长不要超过2048（1024位通常足够）
✅ 理化性质控制在15-20个以内
✅ 如果样本量小于1000，特征总数最好在500以内

如何检测冗余？

有两个方法我经常用：

  1. 相关性热图：计算特征之间的皮尔逊系数。相关系数超过0.85的，只保留其中一个。
  2. 方差选择：方差为0的特征（所有样本值都一样）直接删掉。指纹里有些位可能全是0，毫无信息。
import pandas as pd
from sklearn.feature_selection import VarianceThreshold

# 假设 combined_df 是特征矩阵 (n_samples, n_features)
# 去掉方差为0的特征
selector = VarianceThreshold(threshold=0.0)
combined_reduced = selector.fit_transform(combined_df)

print(f"原始特征数: {combined_df.shape[1]}")
print(f"筛选后特征数: {combined_reduced.shape[1]}")
我曾经在项目早期因为没做这一步，用2060维特征拟合了一个400样本的训练集。结果XGBoost跑了20分钟，训练精度高达0.99，但验证精度只有0.62——典型的过拟合。去掉冗余特征后，特征降到300维，训练时间减少80%，验证精度反而升到0.73。

3.7 小结

分子描述符特征工程不是一次性的工作。它需要你结合业务理解、数据探索和实验验证来反复迭代。

这一章我们聊了：

  • 分子指纹（尤其是ECFP4）的原理和RDKit实现
  • 理化性质的选择与计算
  • 指纹+性质的组合特征
  • 冗余检测与降维

接下来，你就可以把这些特征喂给XGBoost了。但别急——在训练之前，我们还需要聊聊数据分割和评估指标，不然你根本不知道自己的模型好不好。

← 上一章 📖 返回目录 下一章 →