rdkit-cheminformatics

RDKit 化学信息学工具包

概述

RDKit 是 Python 的标准开源化学信息学库，提供全面的 API，用于分子解析、描述符计算、指纹生成、子结构搜索和化学反应。本技能将介绍完整的化合物库分析和虚拟筛选工作流——从加载分子，到类药性过滤、相似性筛选和结果可视化。

适用场景

• 计算一组化合物的分子性质（MW、LogP、TPSA、HBD/HBA）
• 使用指纹相似度筛选库中的化合物，使其与参考化合物进行比较
• 按子结构（SMARTS 模式）筛选化合物，以分析官能团
• 使用 Lipinski 五法则或自定义过滤条件评估类药性
• 生成 2D 示意图或 3D 构象，用于后续对接
• 使用反应 SMARTS 枚举化学库（组合化学）
• 按结构相似性对化合物聚类，以分析多样性
• 标准化并去重分子数据集（规范 SMILES、InChI）
• 若需要支持批处理和错误处理的高阶 RDKit 封装，请改用 datamol-cheminformatics；若需要多格式转换（MOL2、XYZ、PDB），请改用 openbabel

前提条件

• Python 包：rdkit-pypi（或通过 conda 安装的 rdkit）、pandas、matplotlib、numpy
• 数据要求：以 SMILES 字符串、SDF 文件或 MOL 文件表示的分子结构
• 环境：Python 3.8+；完整安装 RDKit 建议使用 conda

# Option 1: pip (lightweight)
pip install rdkit-pypi pandas matplotlib numpy

# Option 2: conda (full features including cartridge)
conda install -c conda-forge rdkit pandas matplotlib numpy

工作流

步骤 1：加载并验证分子

从 SMILES 或 SDF 读取分子结构并验证解析结果。

from rdkit import Chem
import pandas as pd

# --- From SMILES list ---
smiles_list = [
    "CC(=O)Oc1ccccc1C(=O)O",       # Aspirin
    "CC12CCC3C(C1CCC2O)CCC4=CC(=O)CCC34C",  # Testosterone
    "c1ccc2[nH]c(-c3ccccn3)nc2c1",  # Benzimidazole derivative
    "CC(C)Cc1ccc(C(C)C(=O)O)cc1",   # Ibuprofen
    "INVALID_SMILES",                 # Will fail
]

mols = []
failed = []
for smi in smiles_list:
    mol = Chem.MolFromSmiles(smi)
    if mol is not None:
        mol.SetProp("_SMILES", smi)
        mols.append(mol)
    else:
        failed.append(smi)

print(f"Successfully parsed: {len(mols)}/{len(smiles_list)}")
print(f"Failed: {failed}")

# --- From SDF file ---
# suppl = Chem.SDMolSupplier("library.sdf")
# mols = [mol for mol in suppl if mol is not None]
# print(f"Loaded {len(mols)} molecules from SDF")

步骤 2：标准化并去重

将 SMILES 规范化并删除重复项，以确保数据集干净。

from rdkit.Chem.MolStandardize import rdMolStandardize

# Standardize: neutralize charges, remove fragments, canonicalize
uncharger = rdMolStandardize.Uncharger()
chooser = rdMolStandardize.LargestFragmentChooser()

standardized = []
seen_smiles = set()

for mol in mols:
    # Keep largest fragment (remove salts/counterions)
    mol = chooser.choose(mol)
    # Neutralize charges
    mol = uncharger.uncharge(mol)
    # Canonical SMILES for deduplication
    canon_smi = Chem.MolToSmiles(mol)
    if canon_smi not in seen_smiles:
        seen_smiles.add(canon_smi)
        mol.SetProp("canonical_smiles", canon_smi)
        standardized.append(mol)

print(f"After standardization: {len(standardized)} unique molecules")
print(f"Removed {len(mols) - len(standardized)} duplicates/salts")

步骤 3：计算分子描述符

计算每个分子的理化性质。

from rdkit.Chem import Descriptors

records = []
for mol in standardized:
    desc = {
        "SMILES": Chem.MolToSmiles(mol),
        "MW": round(Descriptors.MolWt(mol), 2),
        "LogP": round(Descriptors.MolLogP(mol), 2),
        "TPSA": round(Descriptors.TPSA(mol), 2),
        "HBD": Descriptors.NumHDonors(mol),
        "HBA": Descriptors.NumHAcceptors(mol),
        "RotBonds": Descriptors.NumRotatableBonds(mol),
        "AromaticRings": Descriptors.NumAromaticRings(mol),
        "HeavyAtoms": mol.GetNumHeavyAtoms(),
        "RingCount": Descriptors.RingCount(mol),
    }
    records.append(desc)

df = pd.DataFrame(records)
print(df.to_string(index=False))
print(f"\nDescriptor summary:\n{df.describe().round(2)}")

步骤 4：应用类药性过滤条件

使用 Lipinski 五法则和 Veber 标准筛选化合物。

def lipinski_filter(row):
    """Lipinski Ro5: MW<=500, LogP<=5, HBD<=5, HBA<=10"""
    return (row["MW"] <= 500 and row["LogP"] <= 5 and
            row["HBD"] <= 5 and row["HBA"] <= 10)

def veber_filter(row):
    """Veber: RotBonds<=10, TPSA<=140"""
    return row["RotBonds"] <= 10 and row["TPSA"] <= 140

df["Lipinski"] = df.apply(lipinski_filter, axis=1)
df["Veber"] = df.apply(veber_filter, axis=1)
df["DrugLike"] = df["Lipinski"] & df["Veber"]

print(f"Lipinski pass: {df['Lipinski'].sum()}/{len(df)}")
print(f"Veber pass:    {df['Veber'].sum()}/{len(df)}")
print(f"Drug-like:     {df['DrugLike'].sum()}/{len(df)}")

drug_like_mols = [standardized[i] for i in df[df["DrugLike"]].index]
print(f"\n{len(drug_like_mols)} drug-like compounds retained")

步骤 5：生成指纹并搜索相似性

计算 Morgan 指纹，并以参考化合物为基准进行筛选。

from rdkit.Chem import AllChem
from rdkit import DataStructs

# Reference compound (e.g., known active)
ref_smi = "CC(=O)Oc1ccccc1C(=O)O"  # Aspirin
ref_mol = Chem.MolFromSmiles(ref_smi)
ref_fp = AllChem.GetMorganFingerprintAsBitVect(ref_mol, radius=2, nBits=2048)

# Screen library
results = []
for mol in drug_like_mols:
    fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048)
    tanimoto = DataStructs.TanimotoSimilarity(ref_fp, fp)
    results.append({
        "SMILES": Chem.MolToSmiles(mol),
        "Tanimoto": round(tanimoto, 3),
    })

sim_df = pd.DataFrame(results).sort_values("Tanimoto", ascending=False)
print("Similarity ranking:")
print(sim_df.to_string(index=False))

# Filter by threshold
threshold = 0.3
hits = sim_df[sim_df["Tanimoto"] >= threshold]
print(f"\n{len(hits)} compounds with Tanimoto >= {threshold}")

步骤 6：使用 SMARTS 进行子结构筛选

筛选包含特定官能团的化合物。

# Define SMARTS patterns for functional groups of interest
patterns = {
    "Carboxylic acid": "[CX3](=O)[OX2H1]",
    "Amide": "[CX3](=[OX1])[NX3]",
    "Aromatic ring": "c1ccccc1",
    "Hydroxyl": "[OX2H]",
    "Ester": "[CX3](=O)[OX2][C]",
}

print("Substructure matches:")
for name, smarts in patterns.items():
    query = Chem.MolFromSmarts(smarts)
    match_count = sum(1 for mol in drug_like_mols if mol.HasSubstructMatch(query))
    print(f"  {name}: {match_count}/{len(drug_like_mols)} compounds")

# Get specific matches with atom indices
query = Chem.MolFromSmarts("[CX3](=O)[OX2H1]")  # Carboxylic acid
for mol in drug_like_mols:
    matches = mol.GetSubstructMatches(query)
    if matches:
        smi = Chem.MolToSmiles(mol)
        print(f"\n{smi}: {len(matches)} carboxylic acid group(s)")
        for match in matches:
            print(f"  Atom indices: {match}")

步骤 7：2D 可视化和网格图

生成适用于出版物的分子示意图。

from rdkit.Chem import Draw
from rdkit.Chem.Draw import rdMolDraw2D

# Grid image of top hits
legends = [f"Tan={row['Tanimoto']}" for _, row in sim_df.head(4).iterrows()]
top_mols = [Chem.MolFromSmiles(smi) for smi in sim_df.head(4)["SMILES"]]

img = Draw.MolsToGridImage(
    top_mols,
    molsPerRow=2,
    subImgSize=(300, 300),
    legends=legends,
)
img.save("top_hits_grid.png")
print("Saved top_hits_grid.png")

# Highlight substructure in a molecule
mol = top_mols[0]
query = Chem.MolFromSmarts("[CX3](=O)[OX2H1]")
match = mol.GetSubstructMatch(query)
if match:
    highlight_img = Draw.MolToImage(mol, size=(400, 400), highlightAtoms=match)
    highlight_img.save("substructure_highlight.png")
    print("Saved substructure_highlight.png")

步骤 8：导出结果

保存分析结果和筛选后的化合物。

import os

os.makedirs("results", exist_ok=True)

# Save descriptor table
df.to_csv("results/descriptors.csv", index=False)
print(f"Saved descriptors for {len(df)} compounds to results/descriptors.csv")

# Save drug-like compounds as SDF
writer = Chem.SDWriter("results/drug_like_compounds.sdf")
for i, mol in enumerate(drug_like_mols):
    # Attach descriptors as SDF properties
    row = df[df["DrugLike"]].iloc[i]
    mol.SetProp("MW", str(row["MW"]))
    mol.SetProp("LogP", str(row["LogP"]))
    mol.SetProp("TPSA", str(row["TPSA"]))
    writer.write(mol)
writer.close()
print(f"Saved {len(drug_like_mols)} drug-like compounds to results/drug_like_compounds.sdf")

# Save similarity results
sim_df.to_csv("results/similarity_results.csv", index=False)
print(f"Saved similarity rankings to results/similarity_results.csv")

关键参数

参数	默认值	范围 / 选项	作用
MolFromSmiles(sanitize=)	True	True、False	解析时自动验证并识别芳香性
Morgan radius	2	1-3	指纹半径；2 ≈ ECFP4，3 ≈ ECFP6
Morgan nBits	2048	1024-4096	指纹位数；数值越大，碰撞越少
Tanimoto threshold	0.7	0.3-0.9	相似度阈值；数值越低，筛选越宽松
Lipinski MW cutoff	500	300-600	类药性筛选的最大分子量
Lipinski LogP cutoff	5	3-6	最大亲脂性
Veber RotBonds cutoff	10	7-15	口服生物利用度的最大可旋转键数
Veber TPSA cutoff	140	120-160	最大极性表面积（Ų）
EmbedMolecule(randomSeed=)	None	任意整数	用于可复现 3D 构象生成的种子
Butina distThresh	0.3	0.2-0.5	Butina 聚类的距离阈值

常见用法

用法：使用 Butina 聚类进行多样性筛选

适用场景：从大型化合物库中选出具有多样性的子集。

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

# Generate fingerprints
fps = [AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=2048)
       for mol in standardized]

# Build distance matrix (lower triangle)
dists = []
for i in range(1, len(fps)):
    sims = DataStructs.BulkTanimotoSimilarity(fps[i], fps[:i])
    dists.extend([1 - s for s in sims])

# Cluster
clusters = Butina.ClusterData(dists, len(fps), distThresh=0.3, isDistData=True)
print(f"{len(clusters)} clusters from {len(fps)} compounds")

# Pick centroid from each cluster (first element = centroid)
diverse_indices = [c[0] for c in clusters]
diverse_mols = [standardized[i] for i in diverse_indices]
print(f"Selected {len(diverse_mols)} diverse representatives")

用法：反应枚举（酰胺偶联）

适用场景：通过反应 SMARTS 从构建单元生成组合化学库。

from rdkit.Chem import AllChem, Chem

# Amide coupling: carboxylic acid + amine → amide
rxn = AllChem.ReactionFromSmarts(
    "[C:1](=[O:2])[OH].[N:3]([H])([H])[C:4]>>[C:1](=[O:2])[N:3][C:4]"
)

acids = [Chem.MolFromSmiles(s) for s in ["OC(=O)c1ccccc1", "OC(=O)CC"]]
amines = [Chem.MolFromSmiles(s) for s in ["NCC", "NC1CCCCC1"]]

products = []
for acid in acids:
    for amine in amines:
        ps = rxn.RunReactants((acid, amine))
        for product_set in ps:
            for prod in product_set:
                Chem.SanitizeMol(prod)
                products.append(Chem.MolToSmiles(prod))

print(f"Generated {len(products)} products:")
for p in products:
    print(f"  {p}")

用法：生成 3D 构象并进行 MMFF 优化

适用场景：为分子对接或 3D 药效团分析准备分子。

from rdkit import Chem
from rdkit.Chem import AllChem

mol = Chem.MolFromSmiles("CC(=O)Oc1ccccc1C(=O)O")
mol = Chem.AddHs(mol)  # Required for 3D embedding

# Generate multiple conformers
params = AllChem.ETKDGv3()
params.randomSeed = 42
params.numThreads = 0  # Use all available cores
conf_ids = AllChem.EmbedMultipleConfs(mol, numConfs=10, params=params)
print(f"Generated {len(conf_ids)} conformers")

# Optimize with MMFF94 force field
energies = []
for conf_id in conf_ids:
    result = AllChem.MMFFOptimizeMolecule(mol, confId=conf_id)
    ff = AllChem.MMFFGetMoleculeForceField(mol, AllChem.MMFFGetMoleculeProperties(mol), confId=conf_id)
    energy = ff.CalcEnergy()
    energies.append((conf_id, energy))
    print(f"  Conformer {conf_id}: {energy:.2f} kcal/mol (converged={result == 0})")

# Get lowest energy conformer
best_id = min(energies, key=lambda x: x[1])[0]
print(f"\nBest conformer: {best_id} ({min(e for _, e in energies):.2f} kcal/mol)")

# Save to SDF
writer = Chem.SDWriter("conformers.sdf")
for conf_id, energy in energies:
    mol.SetProp("Energy", f"{energy:.2f}")
    writer.write(mol, confId=conf_id)
writer.close()

配方：带原子索引和自定义绘图的分子可视化

适用场景：调试 SMARTS 匹配，为报告标注原子位置。

from rdkit import Chem
from rdkit.Chem.Draw import rdMolDraw2D

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

# Custom drawer with atom indices and stereo annotations
drawer = rdMolDraw2D.MolDraw2DCairo(500, 400)
opts = drawer.drawOptions()
opts.addAtomIndices = True
opts.addStereoAnnotation = True
opts.bondLineWidth = 2.0

drawer.DrawMolecule(mol)
drawer.FinishDrawing()

with open("annotated_molecule.png", "wb") as f:
    f.write(drawer.GetDrawingText())
print("Saved annotated_molecule.png with atom indices")

预期输出

• results/descriptors.csv — 描述符表格（SMILES、MW、LogP、TPSA、HBD、HBA、RotBonds、Lipinski、Veber、DrugLike）
• results/drug_like_compounds.sdf — 经过筛选且附带属性的 SDF 格式化合物
• results/similarity_results.csv — 与参考化合物进行比较的 Tanimoto 相似度排名
• top_hits_grid.png — 相似度最高的化合物的 2D 网格图像
• substructure_highlight.png — 突出显示官能团的分子图像
• conformers.sdf — 带 MMFF 能量的 3D 构象集合
• annotated_molecule.png — 带原子索引的 2D 示意图

故障排除

问题	原因	解决方案
MolFromSmiles 返回 None	SMILES 字符串无效或存在价态错误	检查 SMILES 是否有效；使用 Chem.MolFromSmiles(smi, sanitize=False)，然后调用 DetectChemistryProblems() 进行诊断
Kekulization error	芳香环系统无效	检查是否存在非标准芳香性；尝试 Chem.SanitizeMol(mol, sanitizeOps=Chem.SANITIZE_ALL ^ Chem.SANITIZE_KEKULIZE)
EmbedMolecule 返回 -1	3D 嵌入失败（环张力、空间位阻）	使用 AllChem.EmbedMolecule(mol, maxAttempts=50, useRandomCoords=True)
MMFF has null force field	缺少原子类型对应的 MMFF 参数	切换到 UFF：AllChem.UFFOptimizeMolecule(mol)
描述符值错误	缺少显式氢	对依赖氢的属性计算描述符前，调用 Chem.AddHs(mol)
ForwardSDMolSupplier 为空	文件不存在或格式错误	验证路径；使用 Chem.SDMolSupplier(path, sanitize=False) 跳过有问题的分子
大型化合物库的指纹计算缓慢	顺序处理	使用预分配数组调用 AllChem.GetMorganFingerprintAsBitVect；考虑使用 rdkit.Chem.MultithreadedSDMolSupplier
SMARTS 模式未匹配到任何内容	SMARTS 语法错误或芳香性不匹配	先用简单的 SMILES 测试模式；表示任意碳原子时使用 [#6] 而非 C
大型 SDF 出现 MemoryError	将整个文件加载到内存	使用 ForwardSDMolSupplier 进行流式读取；分批处理
规范 SMILES 不一致	RDKit 版本不同	固定 RDKit 版本；显式使用 Chem.MolToSmiles(mol, canonical=True)

随附资源

此技能在 references/ 子目录中包含参考文件：

• references/api_reference.md — 按功能分类整理的关键 RDKit 模块与函数索引（I/O、描述符、指纹、绘图、反应）
• references/descriptors_guide.md — 200 多种分子描述符的完整列表，包含名称、说明和典型范围
• references/smarts_patterns.md — 用于检测官能团的常见 SMARTS 模式，按化学类别整理

参考资料

• RDKit 文档 — 官方文档和入门指南
• RDKit Cookbook — 配方和常见工作流模式
• Python 中的 RDKit 入门 — 全面教程
• Lipinski, C.A. 等人（2001）Advanced Drug Delivery Reviews 46:3-26 — 五规则
• Veber, D.F. 等人（2002）J. Med. Chem. 45:2615-2623 — 口服生物利用度标准
• Morgan, H.L.（1965）J. Chem. Doc. 5:107-113 — 环状指纹算法