Python 实战：基于随机森林与 XGBoost 的 729 维药物分子描述符特征筛选

1. 高维特征筛选的挑战与方法论

面对729维分子描述符的庞大规模，传统统计方法往往束手无策。1974个样本与729个特征的比值约为2.7:1，远低于机器学习领域推荐的10:1样本特征比。这种高维稀疏场景下，我们需要一套组合拳式的特征筛选策略：

低方差滤波 作为第一道防线，可快速剔除225个零方差特征。这类似于建筑前的土地平整，为后续精细筛选创造条件。在Python中实现仅需几行代码：

from sklearn.feature_selection import VarianceThreshold
selector = VarianceThreshold(threshold=0)
selector.fit(molecular_descriptors)
remaining_features = selector.get_support(indices=True)

灰色关联分析 （GRA）则像分子世界的"相亲大会"，计算每个特征与pIC50的关联度。分辨系数设置为0.5时，前200名特征的得分差距达0.23，而201名后仅差0.068——这个陡峭的断层暗示着关键生物活性因子的集中区域。

2. 递归特征消除的工程实践

随机森林的递归特征消除(RFE-RF)是本次筛选的核心武器。与传统RFE不同，我们采用50次重复实验的稳健策略：

from sklearn.ensemble import RandomForestRegressor
from sklearn.feature_selection import RFE

rfe_results = []
for _ in range(50):
    estimator = RandomForestRegressor(n_estimators=20, max_depth=4)
    selector = RFE(estimator, n_features_to_select=30, step=1)
    selector.fit(X_gra, y_pIC50)
    rfe_results.append(selector.support_)

这种蒙特卡洛式的实验设计有效抵消了算法随机性，最终得到的30个高频特征构成"精英俱乐部"。进一步通过10次随机森林回归验证，最终锁定20个核心描述符。

关键发现：前20个特征在MIC、距离相关系数和随机森林重要性三个维度上表现高度一致，但Pearson相关系数均低于0.6，证实生物活性与分子特征间存在显著非线性关系。

3. 特征独立性优化算法

Top30特征间可能存在冗余信息。我们设计了一种类NMS（非极大值抑制）算法：

1. 计算30×30距离相关系数矩阵
2. 按特征重要性降序处理
3. 剔除相关系数>0.6的次重要特征
4. 保留特征子集信息量最大化

def nms_feature_selection(features, scores, threshold=0.6):
    # 计算距离相关系数矩阵
    corr_matrix = pairwise_distances(features.T, metric=distance_correlation)
    # 按得分排序
    idx = np.argsort(scores)[::-1]
    selected = []
    
    while len(idx) > 0:
        current = idx[0]
        selected.append(current)
        # 找出强相关特征
        correlated = np.where(corr_matrix[current] > threshold)[0]
        # 移除已选特征及其强相关特征
        idx = np.setdiff1d(idx, np.union1d([current], correlated))
    
    return selected

4. 模型构建与性能对比

在筛选后的特征空间上，我们对比了五种非线性模型：

GBDT以R²=0.807胜出，其超参优化过程展示了典型的学习曲线：

param_grid = {
    'n_estimators': range(240,300),
    'max_depth': range(1,14),
    'subsample': [0.6,0.7,0.8,0.9,1]
}

best_params = {'n_estimators': 295, 'max_depth': 5, 'subsample': 0.7}

模型预测值与真实值的散点图显示，大多数样本落在对角线附近，但在pIC50极高和极低区域存在轻微系统性偏差——这提示我们可能需要引入更复杂的分子相互作用特征。

5. 关键分子描述符解析

最终筛选出的Top20特征中，有三个典型代表值得深入探讨：

1. nRing （环数量）：与pIC50呈正相关（MIC=0.62），暗示多环结构可能增强与ERα受体的结合稳定性
2. SpMax_B（m） （最大自相关指数）：反映分子电子分布特性，距离相关系数达0.58
3. GATS2p （Geary自相关）：描述分子拓扑结构，在50次RFE中出现频率高达47次

这些发现与已知的ERα拮抗剂结构特征高度吻合——多环芳烃结构和特定电子分布确实是药物设计中重点考虑的因素。

6. 工程实现建议

对于实际药物研发项目，建议采用以下优化流程：

1. 数据预处理

   • 分子描述符标准化
   • 活性值(pIC50)归一化
   • 缺失值插补（中位数/模型预测）

2. 特征筛选流水线

from sklearn.pipeline import Pipeline

pipeline = Pipeline([
    ('variance_threshold', VarianceThreshold(threshold=0)),
    ('gra_selection', GRASelector(top_k=200)),
    ('rfe', RFE(estimator=RandomForestRegressor(), n_features_to_select=30)),
    ('final_selection', NMSSelector(threshold=0.6))
])

3. 模型部署注意事项
   • 定期用新化合物数据更新模型
   • 建立特征漂移监测机制
   • 开发SHAP值解释模块增强可解释性

这种结构化特征选择方案不仅适用于药物发现，也可迁移到材料设计、化工产品开发等其他高维特征筛选场景。关键在于根据具体问题调整各阶段的阈值参数，就像调整显微镜的焦距，直到关键特征清晰呈现。