rdkit

RDKit - 化学信息学与药物发现

RDKit 是化学信息学领域的行业标准开源工具包。它提供全面的分子操作、描述符计算、指纹生成、子结构搜索和 3D 分子建模工具。制药公司广泛使用 RDKit 进行药物发现和虚拟筛选。

适用场景

• 读取和写入化学文件格式（SMILES、SDF、MOL2、PDB）。
• 计算分子描述符和类药物性质（Lipinski 五法则）。
• 生成分子指纹以进行相似性搜索。
• 子结构搜索和化学模式匹配（SMARTS）。
• 生成 3D 构象并进行分子对齐。
• 对化合物库进行虚拟筛选。
• 药效团建模和形状相似性分析。
• QSAR（定量构效关系）建模。
• 反应枚举和逆合成。
• 将化学结构可视化为 2D 和 3D 图像。
• 构建用于预测分子性质的机器学习模型。

参考文档

官方文档：https://www.rdkit.org/docs/
RDKit 手册：https://www.rdkit.org/docs/RDKit_Book.html
GitHub：https://github.com/rdkit/rdkit
搜索模式：rdkit.Chem、rdkit.Chem.Descriptors、rdkit.Chem.AllChem、rdkit.DataStructs

核心原则

分子表示

RDKit 将分子表示为图，其中原子是节点，化学键是边。核心对象是 Mol，可通过 SMILES、SDF 文件创建，也可通过编程方式构建。

SMILES（简化分子线性输入规范）

一种基于文本的化学结构表示法。例如：CCO 表示乙醇，c1ccccc1 表示苯。RDKit 可以解析并生成 SMILES 字符串。

用于相似性分析的指纹

分子指纹是编码结构特征的二进制向量。它们可用于对大型化合物库进行快速相似性搜索和聚类。

惰性求值

RDKit 的许多操作采用惰性求值方式——仅在需要时才计算属性。这让针对大型化合物库的操作效率很高。

快速参考

安装

# Via conda (recommended)
conda install -c conda-forge rdkit

# Via pip
pip install rdkit

# For visualization
pip install rdkit pillow

标准导入

from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, Draw, Lipinski
from rdkit.Chem import rdFingerprintGenerator
from rdkit import DataStructs
import numpy as np
import pandas as pd

基本模式 - 从 SMILES 创建分子

from rdkit import Chem

# 1. Create molecule from SMILES
smiles = "CC(=O)OC1=CC=CC=C1C(=O)O"  # Aspirin
mol = Chem.MolFromSmiles(smiles)

# 2. Check if molecule is valid
if mol is None:
    print("Invalid SMILES")
else:
    print(f"Molecular formula: {Chem.rdMolDescriptors.CalcMolFormula(mol)}")
    print(f"Molecular weight: {Descriptors.MolWt(mol):.2f}")

基本模式 - 计算性质

from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski

mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O")

# Calculate drug-like properties
mw = Descriptors.MolWt(mol)
logp = Descriptors.MolLogP(mol)
hbd = Lipinski.NumHDonors(mol)
hba = Lipinski.NumHAcceptors(mol)

print(f"MW: {mw:.2f}, LogP: {logp:.2f}, HBD: {hbd}, HBA: {hba}")

# Check Lipinski's Rule of Five
lipinski_pass = (mw <= 500 and logp <= 5 and hbd <= 5 and hba <= 10)
print(f"Lipinski compliant: {lipinski_pass}")

关键规则

✅ 应该做

• 始终验证分子 - 解析 SMILES/文件后检查 mol is not None，以捕获无效结构。
• 使用规范 SMILES - 使用 Chem.MolToSmiles(mol) 获取规范化（标准化）的 SMILES，以便进行比较。
• 清理分子 - RDKit 默认会自动清理分子（检查价态、芳香性）。除非有特定原因，否则请保持启用。
• 生成 3D 坐标 - 在对齐或对接等 3D 操作前，使用 AllChem.EmbedMolecule()。
• 对大型化合物库使用指纹 - 对数百万种化合物进行相似性搜索时，指纹比直接比较快 1000 倍。
• 指定随机种子 - 为确保构象生成可复现，请始终设置 randomSeed。
• 处理立体化学 - 使用 Chem.AssignStereochemistry() 正确指定 R/S 和 E/Z 标记。
• 批量处理 - 处理数百万个分子时，使用生成器或分块处理，以避免内存问题。

❌ 不应该做

• 不要忽略无效分子 - 始终处理 MolFromSmiles() 返回 None 的情况。
• 不要直接比较 SMILES 字符串 - 不同的 SMILES 可能表示同一个分子。请使用规范 SMILES 或 InChI。
• 不要跳过凯库勒化 - 对于芳香体系，确保正确分配凯库勒结构。
• 不要用描述符进行相似性分析 - 应使用指纹（速度更快，也更合适）。
• 不要忘记氢原子 - 进行 3D 操作时，如有需要，请使用 Chem.AddHs() 添加显式氢原子。
• 不要过度使用 3D 最小化 - 能量最小化速度较慢；仅在必要时使用（如对接、可视化）。

反模式（绝对不要）

from rdkit import Chem
from rdkit.Chem import AllChem

# ❌ BAD: Not checking if molecule is valid
smiles = "INVALID_SMILES"
mol = Chem.MolFromSmiles(smiles)
mw = Descriptors.MolWt(mol)  # Crashes!

# ✅ GOOD: Always validate
mol = Chem.MolFromSmiles(smiles)
if mol is not None:
    mw = Descriptors.MolWt(mol)
else:
    print("Invalid SMILES")

# ❌ BAD: Comparing SMILES strings directly
smiles1 = "CC(C)C"  # isobutane
smiles2 = "C(C)CC"  # same molecule, different SMILES
if smiles1 == smiles2:  # False, but same molecule!
    print("Same")

# ✅ GOOD: Use canonical SMILES
mol1 = Chem.MolFromSmiles(smiles1)
mol2 = Chem.MolFromSmiles(smiles2)
can1 = Chem.MolToSmiles(mol1)
can2 = Chem.MolToSmiles(mol2)
if can1 == can2:  # True
    print("Same molecule")

# ❌ BAD: 3D operations without 3D coordinates
mol = Chem.MolFromSmiles("CCO")
AllChem.AlignMol(mol, ref_mol)  # Fails! No 3D coords

# ✅ GOOD: Generate 3D coordinates first
mol = Chem.MolFromSmiles("CCO")
AllChem.EmbedMolecule(mol)
AllChem.AlignMol(mol, ref_mol)

分子 I/O 与转换

SMILES 解析

from rdkit import Chem

# Parse SMILES
mol = Chem.MolFromSmiles("CCO")

# Parse SMILES with sanitization control
mol = Chem.MolFromSmiles("CCO", sanitize=True)  # Default

# Generate canonical SMILES
canonical = Chem.MolToSmiles(mol)

# Generate isomeric SMILES (includes stereochemistry)
iso_smiles = Chem.MolToSmiles(mol, isomericSmiles=True)

# Generate SMILES without stereochemistry
non_iso = Chem.MolToSmiles(mol, isomericSmiles=False)

# Handle invalid SMILES
smiles_list = ["CCO", "INVALID", "c1ccccc1"]
mols = []
for smi in smiles_list:
    mol = Chem.MolFromSmiles(smi)
    if mol is not None:
        mols.append(mol)
    else:
        print(f"Failed to parse: {smi}")

读取 SDF 文件

from rdkit import Chem

# Read single molecule from file
mol = Chem.MolFromMolFile("molecule.mol")

# Read multiple molecules from SDF
suppl = Chem.SDMolSupplier("compounds.sdf")

# Iterate through molecules
for mol in suppl:
    if mol is None:
        continue
    
    smiles = Chem.MolToSmiles(mol)
    print(f"SMILES: {smiles}")
    
    # Access SDF properties
    if mol.HasProp("_Name"):
        name = mol.GetProp("_Name")
        print(f"Name: {name}")

# Read with removeHs=False to keep explicit hydrogens
suppl = Chem.SDMolSupplier("compounds.sdf", removeHs=False)

写入 SDF 文件

from rdkit import Chem

# Write single molecule
mol = Chem.MolFromSmiles("CCO")
writer = Chem.SDWriter("output.sdf")
writer.write(mol)
writer.close()

# Write multiple molecules
mols = [Chem.MolFromSmiles(s) for s in ["CCO", "c1ccccc1", "CC(=O)O"]]
writer = Chem.SDWriter("output.sdf")
for mol in mols:
    if mol is not None:
        writer.write(mol)
writer.close()

# Add properties to molecules
mol = Chem.MolFromSmiles("CCO")
mol.SetProp("_Name", "Ethanol")
mol.SetProp("Activity", "10.5")
writer = Chem.SDWriter("output.sdf")
writer.write(mol)
writer.close()

InChI 和 InChIKey

from rdkit import Chem

mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O")  # Aspirin

# Generate InChI (unique chemical identifier)
inchi = Chem.MolToInchi(mol)
print(f"InChI: {inchi}")

# Generate InChIKey (hashed InChI, good for database lookups)
inchikey = Chem.MolToInchiKey(mol)
print(f"InChIKey: {inchikey}")

# Parse InChI
mol_from_inchi = Chem.MolFromInchi(inchi)

分子描述符

常用描述符

from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski, Crippen

mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O")  # Aspirin

# Basic properties
mw = Descriptors.MolWt(mol)
num_atoms = mol.GetNumAtoms()
num_heavy_atoms = Lipinski.HeavyAtomCount(mol)

# Lipinski's Rule of Five parameters
logp = Descriptors.MolLogP(mol)
hbd = Lipinski.NumHDonors(mol)
hba = Lipinski.NumHAcceptors(mol)
rotatable_bonds = Lipinski.NumRotatableBonds(mol)

# Topological descriptors
tpsa = Descriptors.TPSA(mol)  # Topological polar surface area
rings = Lipinski.RingCount(mol)
aromatic_rings = Lipinski.NumAromaticRings(mol)

# Complexity
bertz_ct = Descriptors.BertzCT(mol)  # Molecular complexity

print(f"""
Molecular Weight: {mw:.2f}
LogP: {logp:.2f}
HBD: {hbd}
HBA: {hba}
TPSA: {tpsa:.2f}
Rotatable Bonds: {rotatable_bonds}
Aromatic Rings: {aromatic_rings}
""")

计算所有描述符

from rdkit import Chem
from rdkit.Chem import Descriptors

mol = Chem.MolFromSmiles("CCO")

# Get all available descriptors
descriptor_names = [desc[0] for desc in Descriptors.descList]

# Calculate all descriptors
descriptors = {}
for name in descriptor_names:
    calc = getattr(Descriptors, name)
    descriptors[name] = calc(mol)

print(f"Total descriptors: {len(descriptors)}")
print(f"First 5: {list(descriptors.items())[:5]}")

药物相似性筛选

from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski

def check_lipinski(mol):
    """Check Lipinski's Rule of Five."""
    mw = Descriptors.MolWt(mol)
    logp = Descriptors.MolLogP(mol)
    hbd = Lipinski.NumHDonors(mol)
    hba = Lipinski.NumHAcceptors(mol)
    
    rules = {
        'MW <= 500': mw <= 500,
        'LogP <= 5': logp <= 5,
        'HBD <= 5': hbd <= 5,
        'HBA <= 10': hba <= 10
    }
    
    passed = all(rules.values())
    return passed, rules

def check_veber(mol):
    """Check Veber's rules for oral bioavailability."""
    rotatable = Lipinski.NumRotatableBonds(mol)
    tpsa = Descriptors.TPSA(mol)
    
    rules = {
        'Rotatable bonds <= 10': rotatable <= 10,
        'TPSA <= 140': tpsa <= 140
    }
    
    passed = all(rules.values())
    return passed, rules

# Usage
mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O")  # Aspirin
lipinski_pass, lipinski_rules = check_lipinski(mol)
veber_pass, veber_rules = check_veber(mol)

print(f"Lipinski: {lipinski_pass}")
print(f"Veber: {veber_pass}")

分子指纹

Morgan 指纹（圆形）

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import DataStructs

# Create molecule
mol = Chem.MolFromSmiles("CCO")

# Generate Morgan fingerprint (ECFP4)
fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048)

# Convert to numpy array
import numpy as np
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)

# Generate count-based fingerprint (for feature importance)
fp_counts = AllChem.GetMorganFingerprint(mol, radius=2)

# Get feature info
info = {}
fp_info = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048, bitInfo=info)
print(f"Number of on-bits: {len(info)}")

指纹相似度

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import DataStructs

# Create molecules
mol1 = Chem.MolFromSmiles("CCO")
mol2 = Chem.MolFromSmiles("CCCO")
mol3 = Chem.MolFromSmiles("c1ccccc1")

# Generate fingerprints
fp1 = AllChem.GetMorganFingerprintAsBitVect(mol1, radius=2)
fp2 = AllChem.GetMorganFingerprintAsBitVect(mol2, radius=2)
fp3 = AllChem.GetMorganFingerprintAsBitVect(mol3, radius=2)

# Calculate Tanimoto similarity
sim_12 = DataStructs.TanimotoSimilarity(fp1, fp2)
sim_13 = DataStructs.TanimotoSimilarity(fp1, fp3)

print(f"Ethanol vs Propanol: {sim_12:.3f}")  # High similarity
print(f"Ethanol vs Benzene: {sim_13:.3f}")   # Low similarity

# Calculate Dice similarity
dice_12 = DataStructs.DiceSimilarity(fp1, fp2)

# Bulk similarity (compare one to many)
fps = [fp1, fp2, fp3]
similarities = DataStructs.BulkTanimotoSimilarity(fp1, fps)
print(f"Bulk similarities: {similarities}")

其他指纹类型

from rdkit import Chem
from rdkit.Chem import AllChem, RDKFingerprint

mol = Chem.MolFromSmiles("CCO")

# RDKit fingerprint (topological)
fp_rdkit = Chem.RDKFingerprint(mol)

# Atom pair fingerprint
fp_atompair = AllChem.GetHashedAtomPairFingerprintAsBitVect(mol)

# Topological torsion fingerprint
fp_torsion = AllChem.GetHashedTopologicalTorsionFingerprintAsBitVect(mol)

# MACCS keys (166-bit structural keys)
from rdkit.Chem import MACCSkeys
fp_maccs = MACCSkeys.GenMACCSKeys(mol)

子结构搜索

SMARTS 模式匹配

from rdkit import Chem

# Define molecule
mol = Chem.MolFromSmiles("CC(=O)Oc1ccccc1C(=O)O")  # Aspirin

# SMARTS pattern for carboxylic acid
pattern = Chem.MolFromSmarts("C(=O)O")

# Check if substructure exists
has_match = mol.HasSubstructMatch(pattern)
print(f"Contains carboxylic acid: {has_match}")

# Get matching atoms
matches = mol.GetSubstructMatches(pattern)
print(f"Number of matches: {len(matches)}")
print(f"Matching atom indices: {matches}")

# Common SMARTS patterns
patterns = {
    'Carboxylic acid': 'C(=O)O',
    'Ester': 'C(=O)O[C,c]',
    'Amide': 'C(=O)N',
    'Alcohol': '[OH][C,c]',
    'Primary amine': '[NH2][C,c]',
    'Aromatic ring': 'c1ccccc1'
}

for name, smarts in patterns.items():
    pattern = Chem.MolFromSmarts(smarts)
    if mol.HasSubstructMatch(pattern):
        print(f"Contains {name}")

子结构筛选

from rdkit import Chem

# Library of molecules
smiles_list = [
    "CC(=O)O",           # Acetic acid
    "CCO",               # Ethanol
    "c1ccccc1C(=O)O",    # Benzoic acid
    "CCCC",              # Butane
]

mols = [Chem.MolFromSmiles(s) for s in smiles_list]

# Filter molecules containing carboxylic acid
pattern = Chem.MolFromSmarts("C(=O)O")
filtered = [mol for mol in mols if mol.HasSubstructMatch(pattern)]

print(f"Molecules with carboxylic acid: {len(filtered)}/{len(mols)}")

# Filter by multiple patterns (AND logic)
pattern1 = Chem.MolFromSmarts("c1ccccc1")  # Aromatic ring
pattern2 = Chem.MolFromSmarts("C(=O)O")    # Carboxylic acid

aromatic_acids = [
    mol for mol in mols
    if mol.HasSubstructMatch(pattern1) and mol.HasSubstructMatch(pattern2)
]

替换子结构

from rdkit import Chem
from rdkit.Chem import AllChem

# Replace carboxylic acid with ester
mol = Chem.MolFromSmiles("CC(=O)O")  # Acetic acid

# Define replacement
rxn = AllChem.ReactionFromSmarts('[C:1](=O)O>>[C:1](=O)OC')

# Apply reaction
products = rxn.RunReactants((mol,))

if products:
    product = products[0][0]
    print(f"Product: {Chem.MolToSmiles(product)}")  # Methyl acetate

3D 构象生成

生成 3D 坐标

from rdkit import Chem
from rdkit.Chem import AllChem

# Create molecule
mol = Chem.MolFromSmiles("CCO")

# Add hydrogens (required for 3D)
mol = Chem.AddHs(mol)

# Generate 3D coordinates
result = AllChem.EmbedMolecule(mol, randomSeed=42)

if result == 0:  # Success
    print("3D coordinates generated")
else:
    print("Failed to generate 3D coordinates")

# Optimize geometry with MMFF force field
AllChem.MMFFOptimizeMolecule(mol)

# Get atomic positions
conf = mol.GetConformer()
for i in range(mol.GetNumAtoms()):
    pos = conf.GetAtomPosition(i)
    print(f"Atom {i}: ({pos.x:.3f}, {pos.y:.3f}, {pos.z:.3f})")

多个构象

from rdkit import Chem
from rdkit.Chem import AllChem

mol = Chem.MolFromSmiles("CCCC")  # Butane
mol = Chem.AddHs(mol)

# Generate multiple conformers
conf_ids = AllChem.EmbedMultipleConfs(
    mol,
    numConfs=10,
    randomSeed=42,
    pruneRmsThresh=0.5  # Remove similar conformers
)

print(f"Generated {len(conf_ids)} conformers")

# Optimize each conformer
for conf_id in conf_ids:
    AllChem.MMFFOptimizeMolecule(mol, confId=conf_id)

# Get energies
props = AllChem.MMFFGetMoleculeProperties(mol)
for conf_id in conf_ids:
    ff = AllChem.MMFFGetMoleculeForceField(mol, props, confId=conf_id)
    energy = ff.CalcEnergy()
    print(f"Conformer {conf_id}: {energy:.2f} kcal/mol")

分子对齐

from rdkit import Chem
from rdkit.Chem import AllChem

# Reference molecule
ref_mol = Chem.MolFromSmiles("c1ccccc1C")  # Toluene
ref_mol = Chem.AddHs(ref_mol)
AllChem.EmbedMolecule(ref_mol)

# Probe molecule
probe_mol = Chem.MolFromSmiles("c1ccccc1CC")  # Ethylbenzene
probe_mol = Chem.AddHs(probe_mol)
AllChem.EmbedMolecule(probe_mol)

# Align probe to reference
rmsd = AllChem.AlignMol(probe_mol, ref_mol)
print(f"RMSD: {rmsd:.3f} Å")

# Get aligned coordinates
# Now probe_mol has coordinates aligned to ref_mol

分子可视化

2D 绘图

from rdkit import Chem
from rdkit.Chem import Draw
import matplotlib.pyplot as plt

# Single molecule
mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O")  # Aspirin
img = Draw.MolToImage(mol, size=(300, 300))
plt.imshow(img)
plt.axis('off')
plt.show()

# Multiple molecules
mols = [Chem.MolFromSmiles(s) for s in ["CCO", "c1ccccc1", "CC(=O)O"]]
legends = ["Ethanol", "Benzene", "Acetic acid"]

img = Draw.MolsToGridImage(
    mols,
    molsPerRow=3,
    subImgSize=(200, 200),
    legends=legends
)
plt.imshow(img)
plt.axis('off')
plt.show()

高亮子结构

from rdkit import Chem
from rdkit.Chem import Draw

mol = Chem.MolFromSmiles("CC(=O)OC1=CC=CC=C1C(=O)O")  # Aspirin

# Highlight carboxylic acid group
pattern = Chem.MolFromSmarts("C(=O)O")
match = mol.GetSubstructMatch(pattern)

# Draw with highlighted atoms
img = Draw.MolToImage(mol, highlightAtoms=match, size=(300, 300))

保存到文件

from rdkit import Chem
from rdkit.Chem import Draw

mol = Chem.MolFromSmiles("CCO")

# Save as PNG
Draw.MolToFile(mol, "molecule.png", size=(300, 300))

# Save as SVG (vector graphics)
from rdkit.Chem.Draw import rdMolDraw2D

drawer = rdMolDraw2D.MolDraw2DSVG(300, 300)
drawer.DrawMolecule(mol)
drawer.FinishDrawing()
svg = drawer.GetDrawingText()

with open("molecule.svg", "w") as f:
    f.write(svg)

实用工作流

1. 虚拟筛选流程

from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors, Lipinski
from rdkit import DataStructs
import pandas as pd

def screen_library(library_file, reference_smiles, similarity_threshold=0.7):
    """Screen compound library for similar, drug-like molecules."""
    
    # Reference molecule and fingerprint
    ref_mol = Chem.MolFromSmiles(reference_smiles)
    ref_fp = AllChem.GetMorganFingerprintAsBitVect(ref_mol, radius=2)
    
    # Results
    hits = []
    
    # Read library
    suppl = Chem.SDMolSupplier(library_file)
    
    for i, mol in enumerate(suppl):
        if mol is None:
            continue
        
        # Step 1: Drug-likeness filter (Lipinski)
        mw = Descriptors.MolWt(mol)
        logp = Descriptors.MolLogP(mol)
        hbd = Lipinski.NumHDonors(mol)
        hba = Lipinski.NumHAcceptors(mol)
        
        if not (mw <= 500 and logp <= 5 and hbd <= 5 and hba <= 10):
            continue
        
        # Step 2: Similarity filter
        fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2)
        similarity = DataStructs.TanimotoSimilarity(ref_fp, fp)
        
        if similarity < similarity_threshold:
            continue
        
        # Step 3: PAINS filter (Pan-Assay Interference Compounds)
        from rdkit.Chem import FilterCatalog
        params = FilterCatalog.FilterCatalogParams()
        params.AddCatalog(FilterCatalog.FilterCatalogParams.FilterCatalogs.PAINS)
        catalog = FilterCatalog.FilterCatalog(params)
        
        if catalog.HasMatch(mol):
            continue
        
        # Passed all filters
        hits.append({
            'id': i,
            'smiles': Chem.MolToSmiles(mol),
            'similarity': similarity,
            'mw': mw,
            'logp': logp
        })
    
    # Convert to DataFrame
    df_hits = pd.DataFrame(hits)
    df_hits = df_hits.sort_values('similarity', ascending=False)
    
    return df_hits

# Usage
# hits = screen_library('compounds.sdf', 'CC(=O)OC1=CC=CC=C1C(=O)O')
# print(f"Found {len(hits)} hits")

2. 多样性选择

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

def select_diverse_set(smiles_list, n_select=100):
    """Select diverse subset using MaxMin algorithm."""
    
    # Generate fingerprints
    mols = [Chem.MolFromSmiles(s) for s in smiles_list]
    fps = [AllChem.GetMorganFingerprintAsBitVect(m, radius=2) for m in mols if m]
    
    if len(fps) < n_select:
        return list(range(len(fps)))
    
    # MaxMin diversity picking
    from rdkit.SimDivFilters import MaxMinPicker
    
    def distance_function(i, j):
        return 1 - DataStructs.TanimotoSimilarity(fps[i], fps[j])
    
    picker = MaxMinPicker()
    picks = picker.LazyPick(
        distance_function,
        len(fps),
        n_select,
        seed=42
    )
    
    return list(picks)

# Usage
smiles_list = ["CCO", "CCCO", "c1ccccc1", "CC(=O)O", "CCCCCCCC"]
diverse_indices = select_diverse_set(smiles_list, n_select=3)
diverse_smiles = [smiles_list[i] for i in diverse_indices]
print(f"Selected: {diverse_smiles}")

3. QSAR 模型构建

from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
import numpy as np

def build_qsar_model(smiles_list, activities):
    """Build QSAR model from SMILES and activities."""
    
    # Generate fingerprints
    fps = []
    valid_activities = []
    
    for smi, act in zip(smiles_list, activities):
        mol = Chem.MolFromSmiles(smi)
        if mol is not None:
            # Morgan fingerprint as features
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048)
            arr = np.zeros((1,))
            DataStructs.ConvertToNumpyArray(fp, arr)
            fps.append(arr)
            valid_activities.append(act)
    
    X = np.array(fps)
    y = np.array(valid_activities)
    
    # Split data
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=42
    )
    
    # Train model
    model = RandomForestRegressor(n_estimators=100, random_state=42)
    model.fit(X_train, y_train)
    
    # Evaluate
    y_pred_train = model.predict(X_train)
    y_pred_test = model.predict(X_test)
    
    print(f"Train R²: {r2_score(y_train, y_pred_train):.3f}")
    print(f"Test R²: {r2_score(y_test, y_pred_test):.3f}")
    print(f"Test RMSE: {mean_squared_error(y_test, y_pred_test, squared=False):.3f}")
    
    return model

# Usage
# smiles = ["CCO", "CCCO", "c1ccccc1", "CC(=O)O"]
# activities = [5.2, 4.8, 6.1, 5.5]  # pIC50 values
# model = build_qsar_model(smiles, activities)

4. 骨架分析

from rdkit import Chem
from rdkit.Chem.Scaffolds import MurckoScaffold
from collections import Counter

def analyze_scaffolds(smiles_list):
    """Analyze Murcko scaffolds in compound set."""
    
    scaffolds = []
    
    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol is not None:
            # Get Murcko scaffold
            scaffold = MurckoScaffold.GetScaffoldForMol(mol)
            scaffold_smi = Chem.MolToSmiles(scaffold)
            scaffolds.append(scaffold_smi)
    
    # Count scaffolds
    scaffold_counts = Counter(scaffolds)
    
    print(f"Unique scaffolds: {len(scaffold_counts)}")
    print("\nTop 5 scaffolds:")
    for scaffold, count in scaffold_counts.most_common(5):
        print(f"{scaffold}: {count}")
    
    return scaffold_counts

# Usage
# smiles_list = ["c1ccccc1CC", "c1ccccc1CCC", "c1ccc(O)cc1"]
# scaffolds = analyze_scaffolds(smiles_list)

5. 反应枚举

from rdkit import Chem
from rdkit.Chem import AllChem

def enumerate_amide_coupling(acids, amines):
    """Enumerate all possible amide products."""
    
    # Define reaction SMARTS
    rxn = AllChem.ReactionFromSmarts('[C:1](=[O:2])O.[N:3]>>[C:1](=[O:2])[N:3]')
    
    products = []
    
    for acid_smi in acids:
        for amine_smi in amines:
            acid = Chem.MolFromSmiles(acid_smi)
            amine = Chem.MolFromSmiles(amine_smi)
            
            if acid is None or amine is None:
                continue
            
            # Run reaction
            products_tuple = rxn.RunReactants((acid, amine))
            
            if products_tuple:
                product = products_tuple[0][0]
                Chem.SanitizeMol(product)
                product_smi = Chem.MolToSmiles(product)
                products.append({
                    'acid': acid_smi,
                    'amine': amine_smi,
                    'product': product_smi
                })
    
    return products

# Usage
acids = ["CC(=O)O", "c1ccccc1C(=O)O"]
amines = ["CCN", "c1ccccc1N"]
products = enumerate_amide_coupling(acids, amines)
print(f"Generated {len(products)} amides")

性能优化

批量操作

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

# Instead of loop
# ❌ SLOW
# fps = []
# for smi in smiles_list:
#     mol = Chem.MolFromSmiles(smi)
#     fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
#     fps.append(fp)

# ✅ FAST: Use PandasTools for bulk operations
from rdkit.Chem import PandasTools

df = pd.DataFrame({'SMILES': smiles_list})
PandasTools.AddMoleculeColumnToFrame(df, 'SMILES', 'Molecule')

# Calculate properties in bulk
df['MW'] = df['Molecule'].apply(lambda x: Descriptors.MolWt(x) if x else None)
df['LogP'] = df['Molecule'].apply(lambda x: Descriptors.MolLogP(x) if x else None)

并行处理

from rdkit import Chem
from rdkit.Chem import AllChem
from multiprocessing import Pool
import pandas as pd

def process_molecule(smiles):
    """Process single molecule."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    
    return {
        'smiles': smiles,
        'mw': Descriptors.MolWt(mol),
        'logp': Descriptors.MolLogP(mol),
        'fp': AllChem.GetMorganFingerprintAsBitVect(mol, 2)
    }

def process_library_parallel(smiles_list, n_jobs=4):
    """Process library in parallel."""
    with Pool(n_jobs) as pool:
        results = pool.map(process_molecule, smiles_list)
    
    # Filter None results
    results = [r for r in results if r is not None]
    return pd.DataFrame(results)

# Usage
# df = process_library_parallel(large_smiles_list, n_jobs=8)

缓存计算结果

from functools import lru_cache
from rdkit import Chem
from rdkit.Chem import Descriptors

@lru_cache(maxsize=10000)
def get_mol_properties(smiles):
    """Calculate properties with caching."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    
    return {
        'mw': Descriptors.MolWt(mol),
        'logp': Descriptors.MolLogP(mol),
        'tpsa': Descriptors.TPSA(mol)
    }

# Repeated calls will use cache
props1 = get_mol_properties("CCO")  # Calculated
props2 = get_mol_properties("CCO")  # Cached (instant)

常见问题与解决方案

“无效 SMILES”问题

并非所有 SMILES 字符串都有效。

# ❌ Problem: Assuming all SMILES are valid
smiles_list = ["CCO", "INVALID", "c1ccccc1"]
mols = [Chem.MolFromSmiles(s) for s in smiles_list]
# Contains None!

# ✅ Solution: Filter invalid molecules
mols = [Chem.MolFromSmiles(s) for s in smiles_list]
valid_mols = [m for m in mols if m is not None]

# ✅ Better: Track which failed
results = []
for smi in smiles_list:
    mol = Chem.MolFromSmiles(smi)
    if mol is not None:
        results.append({'smiles': smi, 'mol': mol, 'valid': True})
    else:
        results.append({'smiles': smi, 'mol': None, 'valid': False})

“立体化学信息丢失”问题

如果不小心，生成 SMILES 时可能会丢失立体化学信息。

from rdkit import Chem

# Molecule with stereochemistry
chiral_smiles = "C[C@H](O)CC"  # (S)-2-butanol
mol = Chem.MolFromSmiles(chiral_smiles)

# ❌ BAD: Lose stereochemistry
non_iso = Chem.MolToSmiles(mol, isomericSmiles=False)
print(non_iso)  # "CC(O)CC" - lost chirality!

# ✅ GOOD: Preserve stereochemistry
iso = Chem.MolToSmiles(mol, isomericSmiles=True)
print(iso)  # "C[C@H](O)CC" - preserved!

“无氢原子的 3D”问题

3D 操作需要显式氢原子。

from rdkit import Chem
from rdkit.Chem import AllChem

mol = Chem.MolFromSmiles("CCO")

# ❌ BAD: Generate 3D without hydrogens
result = AllChem.EmbedMolecule(mol)
# Poor quality or failure

# ✅ GOOD: Add hydrogens first
mol_h = Chem.AddHs(mol)
result = AllChem.EmbedMolecule(mol_h)
AllChem.MMFFOptimizeMolecule(mol_h)

“指纹类型不匹配”问题

比较不同类型的指纹会得出毫无意义的结果。

from rdkit import Chem
from rdkit.Chem import AllChem, MACCSkeys
from rdkit import DataStructs

mol1 = Chem.MolFromSmiles("CCO")
mol2 = Chem.MolFromSmiles("CCCO")

# ❌ BAD: Comparing different fingerprint types
fp1_morgan = AllChem.GetMorganFingerprintAsBitVect(mol1, 2)
fp2_maccs = MACCSkeys.GenMACCSKeys(mol2)

# This will error or give nonsense!
# similarity = DataStructs.TanimotoSimilarity(fp1_morgan, fp2_maccs)

# ✅ GOOD: Use same fingerprint type
fp1 = AllChem.GetMorganFingerprintAsBitVect(mol1, 2)
fp2 = AllChem.GetMorganFingerprintAsBitVect(mol2, 2)
similarity = DataStructs.TanimotoSimilarity(fp1, fp2)

“内存激增”问题

处理数百万个分子可能会耗尽内存。

# ❌ BAD: Load entire library into memory
suppl = Chem.SDMolSupplier('huge_library.sdf')
mols = [mol for mol in suppl]  # Out of memory!

# ✅ GOOD: Process in batches
def process_in_batches(sdf_file, batch_size=10000):
    suppl = Chem.SDMolSupplier(sdf_file)
    batch = []
    
    for mol in suppl:
        if mol is not None:
            batch.append(mol)
        
        if len(batch) >= batch_size:
            # Process batch
            yield batch
            batch = []
    
    # Process remaining
    if batch:
        yield batch

# Usage
for batch in process_in_batches('huge_library.sdf'):
    # Process each batch
    pass

RDKit 是计算药物发现和化学信息学的基石。它拥有全面的分子操作、描述符计算和相似性搜索工具集，因此对于制药研究、虚拟筛选和化学数据分析不可或缺。掌握 RDKit，你就能通过计算探索广阔的化学空间，并加速药物发现。