rdkit

RDKit - ケモインフォマティクスと創薬

RDKit は、ケモインフォマティクス向けの業界標準のオープンソースツールキットです。分子操作、記述子の計算、フィンガープリントの生成、部分構造検索、3D 分子モデリングのための包括的なツールを提供します。RDKit は、製薬企業で創薬やバーチャルスクリーニングに幅広く利用されています。

使用する場面

• 化学ファイル形式（SMILES、SDF、MOL2、PDB）の読み書き。
• 分子記述子と薬物様特性（リピンスキーのルール・オブ・ファイブ）の計算。
• 類似性検索向けの分子フィンガープリントの生成。
• 部分構造検索と化学パターンマッチング（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 を使用する - 比較用に正準（標準化された）SMILES を取得するには、Chem.MolToSmiles(mol) を使用します。
• 分子をサニタイズする - RDKit はデフォルトで自動的にサニタイズ（原子価チェック、芳香族性の処理）を行います。特別な理由がない限り、有効にしておきます。
• 3D 座標を生成する - アラインメントやドッキングなどの 3D 操作の前に、AllChem.EmbedMolecule() を使用します。
• 大規模なライブラリにはフィンガープリントを使用する - 数百万の化合物の類似性検索では、フィンガープリントを使うと直接比較するより 1000 倍高速です。
• 乱数シードを指定する - 配座生成を再現可能にするため、必ず randomSeed を設定します。
• 立体化学を処理する - R/S および E/Z のラベルを適切に割り当てるには、Chem.AssignStereochemistry() を使用します。
• バッチ処理を行う - メモリの問題を避けるため、数百万の分子を処理する際はジェネレーターやチャンク分割を使用します。

❌ 禁止事項

• 無効な分子を無視しない - MolFromSmiles() が None を返す場合は、必ず適切に処理します。
• SMILES 文字列を直接比較しない - 異なる 2 つの 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)

分子の入出力と変換

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 を習得すれば、広大な化学空間を計算によって探索し、創薬を加速する力を得られます。