分子相似性聚类：如何整理大规模化合物库

聚类把相似的分子归到一起，让「十万个分子」变成「五百个化学系列」。它是库设计、结果去冗余、多样性采样的基础工具——而参数的选择直接决定了聚类是否有意义。

常用的聚类方法

Butina 聚类

from rdkit import Chem, DataStructs
from rdkit.Chem import rdFingerprintGenerator
from rdkit.ML.Cluster import Butina

gen = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=2048)

def butina_cluster(smiles_list, cutoff=0.35):
    """cutoff 是【距离】阈值 = 1 - Tanimoto"""
    mols, valid_smiles = [], []
    for s in smiles_list:
        m = Chem.MolFromSmiles(s)
        if m is not None:
            mols.append(m)
            valid_smiles.append(s)
    fps = [gen.GetFingerprint(m) for m in mols]

    # 【计算下三角距离矩阵】
    dists = []
    for i in range(1, len(fps)):
        sims = DataStructs.BulkTanimotoSimilarity(fps[i], fps[:i])
        dists.extend([1 - s for s in sims])

    clusters = Butina.ClusterData(dists, len(fps), cutoff,
                                  isDistData=True)
    # 【返回的每个簇：第一个元素是「簇中心」】
    print(f"{len(fps)} 个分子 → {len(clusters)} 个簇")
    print(f"最大簇: {len(clusters[0])} 个分子")
    print(f"单例簇: {sum(1 for c in clusters if len(c) == 1)} 个")
    return clusters, valid_smiles

# 【算法原理】：
#   1) 计算每个分子的「邻居数」（距离小于阈值的）
#   2) 邻居最多的作为第一个簇中心
#   3) 它的所有邻居归入该簇，并从候选中移除
#   4) 重复，直到所有分子被分配
#
# 【特点】：
#   - 【簇中心是真实存在的分子】（不是虚拟质心）
#     → 可以直接用作代表
#   - 不需要预先指定簇数
#   - 【贪心算法，结果依赖处理顺序】
#   - 【单例簇（只有 1 个分子）可能很多】

# 【cutoff 的含义】：
#   cutoff = 0.35 表示 Tanimoto > 0.65 的分子聚在一起
#
# 【常用值】：
#   0.2  很严（只有非常相似的才聚在一起）→ 簇多
#   0.35 常用
#   0.4  较松
#   0.5  很松 → 簇少但内部差异大

选择聚类参数

# 【不要凭感觉选 cutoff，看结果决定】

import numpy as np

def scan_cutoff(smiles_list, cutoffs=(0.2, 0.3, 0.35, 0.4, 0.5)):
    """扫描不同 cutoff 的聚类结果"""
    mols = [Chem.MolFromSmiles(s) for s in smiles_list]
    mols = [m for m in mols if m is not None]
    fps = [gen.GetFingerprint(m) for m in mols]
    dists = []
    for i in range(1, len(fps)):
        sims = DataStructs.BulkTanimotoSimilarity(fps[i], fps[:i])
        dists.extend([1 - s for s in sims])

    for c in cutoffs:
        clusters = Butina.ClusterData(dists, len(fps), c, isDistData=True)
        sizes = [len(cl) for cl in clusters]
        singletons = sum(1 for s in sizes if s == 1)
        print(f"cutoff={c:.2f}: {len(clusters):5d} 簇, "
              f"最大 {max(sizes):4d}, "
              f"单例 {singletons:5d} ({singletons/len(clusters):.1%}), "
              f"中位大小 {int(np.median(sizes))}")

# 【选择的判据】：
#   1) 【单例比例】
#      过高（> 70%）→ cutoff 太严，聚类没起作用
#   2) 【最大簇的大小】
#      过大（> 20% 的分子）→ cutoff 太松
#   3) 【簇数与你的需求】
#      如果要选 500 个分子做实验，
#      那么簇数应该在 500 左右
#   4) 【人工检查几个簇】
#      簇内的分子看起来真的相似吗？
#
# 【最实用的做法】：
#   从「我要多少个代表」倒推 cutoff

典型应用

# ---- 应用一：虚拟筛选结果的去冗余 ----
#   对接给出前 5000 个分子，但很多是类似物
#   → 聚类后每簇选最好的

def diverse_selection(smiles_list, scores, n_select=200, cutoff=0.4):
    """聚类后每簇选分数最好的"""
    clusters, valid = butina_cluster(smiles_list, cutoff)
    selected = []
    # 按簇内最好分数排序簇
    cluster_best = []
    for cl in clusters:
        best_idx = max(cl, key=lambda i: scores[i])
        cluster_best.append((best_idx, scores[best_idx]))
    cluster_best.sort(key=lambda x: -x[1])
    for idx, _ in cluster_best[:n_select]:
        selected.append(valid[idx])
    return selected

# 【为什么重要】：
#   如果不去冗余，采购的 200 个化合物
#   可能只覆盖 20 个化学系列
#   → 【实验的信息量大幅降低】

# ---- 应用二：数据划分 ----
#   比骨架划分更严格（见 051）
def cluster_split(smiles_list, frac_train=0.8, cutoff=0.6):
    clusters, valid = butina_cluster(smiles_list, cutoff)
    clusters = sorted(clusters, key=len, reverse=True)
    n_train = frac_train * len(valid)
    train, test = [], []
    for cl in clusters:
        (train if len(train) + len(cl) <= n_train else test).extend(cl)
    return train, test

# ---- 应用三：库的多样性评估 ----
#   簇数 / 分子数 = 多样性指标（见 053）

# ---- 应用四：MaxMin 多样性选择 ----
from rdkit.SimDivFilters import rdSimDivPickers

def maxmin_pick(smiles_list, n_pick=100, seed=42):
    """贪心选择最多样的子集"""
    mols = [Chem.MolFromSmiles(s) for s in smiles_list]
    mols = [m for m in mols if m is not None]
    fps = [gen.GetFingerprint(m) for m in mols]

    picker = rdSimDivPickers.MaxMinPicker()

    def dist_func(i, j):
        return 1 - DataStructs.TanimotoSimilarity(fps[i], fps[j])

    indices = picker.LazyPick(dist_func, len(fps), n_pick, seed=seed)
    return [smiles_list[i] for i in indices]

# 【MaxMin vs 聚类】：
#   MaxMin: 【直接选出 N 个最不相似的】
#     → 适合「我就是要 N 个多样的分子」
#   聚类: 先分组再选代表
#     → 【能看到库的结构（有几个系列、各多大）】
#   → 【想了解库的结构用聚类，只想选样本用 MaxMin】

大规模数据的处理

# 【问题】：Butina 需要完整的距离矩阵
#   n 个分子 → n(n-1)/2 个距离
#   10 万分子 → 50 亿个距离 → 【内存爆炸】
#
# 【应对方案】：
#
# 1) 【分层聚类】
#    先用粗糙但快的方法分成大组
#    （如按骨架、按分子量区间）
#    再在每组内做精细聚类
#
# 2) 【采样后聚类】
#    随机采样 1 万个做聚类，
#    确定簇中心后，把其余分子分配到最近的中心
#
# 3) 【用稀疏的邻居图】
#    只计算相似度高于阈值的对
#    → 用倒排索引或 LSH 加速
#
# 4) 【用专用工具】
#    FPSim2、ChemFP 等做了底层优化
#
# 5) 【降维后用标准聚类】
#    先 PCA/UMAP 降到低维（见 052）
#    再用 K-means 或 HDBSCAN
#    → 【但降维会损失信息，聚类结果的化学意义可能下降】

def hierarchical_cluster_large(smiles_list, first_pass_cutoff=0.6,
                               second_pass_cutoff=0.35, max_group=5000):
    """两阶段聚类处理大数据"""
    # 第一阶段：按骨架粗分
    from rdkit.Chem.Scaffolds import MurckoScaffold
    from collections import defaultdict
    groups = defaultdict(list)
    for i, s in enumerate(smiles_list):
        m = Chem.MolFromSmiles(s)
        if m is None:
            continue
        scaf = MurckoScaffold.MurckoScaffoldSmiles(mol=m) or "acyclic"
        groups[scaf].append(i)

    # 第二阶段：组内精细聚类
    all_clusters = []
    for scaf, indices in groups.items():
        if len(indices) == 1:
            all_clusters.append(indices)
            continue
        if len(indices) > max_group:
            # 太大的组再分
            indices = indices[:max_group]
        sub_smiles = [smiles_list[i] for i in indices]
        clusters, _ = butina_cluster(sub_smiles, second_pass_cutoff)
        for cl in clusters:
            all_clusters.append([indices[k] for k in cl])
    return all_clusters

延伸资源

• Tanimoto：039《Tanimoto 相似性》；骨架：041《Bemis–Murcko Scaffold》；Scaffold Split：051《Scaffold Split》；
• 化学空间可视化：052《化学空间可视化》；多样性评估：053《分子多样性评估》。