rdkit-cheminformatics

RDKit ケモインフォマティクスツールキット

概要

RDKit は Python 用の標準的なオープンソースのケモインフォマティクスライブラリで、分子のパース、記述子の計算、フィンガープリントの生成、部分構造検索、化学反応のための包括的な API を提供します。このスキルでは、化合物ライブラリのプロファイリングとバーチャルスクリーニングのワークフローを、分子の読み込みから薬物様性フィルタリング、類似性スクリーニング、結果の可視化まで、一通り説明します。

使用する場面

• 化合物セットの分子特性（MW、LogP、TPSA、HBD/HBA）を計算する
• フィンガープリントの類似度を使って、参照化合物に対するライブラリのスクリーニングを行う
• 部分構造（SMARTS パターン）で化合物をフィルタリングし、官能基を分析する
• Lipinski の Rule of Five またはカスタムフィルターを使って薬物様性を評価する
• 後続のドッキングに向けて 2D 描画または 3D 配座を生成する
• 反応 SMARTS（コンビナトリアルケミストリー）を使って化学ライブラリを列挙する
• 構造類似性で化合物をクラスタリングし、多様性を分析する
• 分子データセットを標準化し、重複を除去する（canonical 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 の Rule of Five と 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 — 参照化合物に対するタニモト類似度ランキング
• 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 でパターンをテストします。任意の炭素には C ではなく [#6] を使います
大きな SDF で MemoryError が発生する	ファイル全体をメモリに読み込んでいる	ストリーミングには ForwardSDMolSupplier を使い、バッチで処理します
正準 SMILES に一貫性がない	RDKit のバージョンが異なる	RDKit のバージョンを固定します。Chem.MolToSmiles(mol, canonical=True) を明示的に使います

同梱リソース

このスキルには、references/ サブディレクトリに参照ファイルが含まれています。

• references/api_reference.md — 機能別（I/O、記述子、フィンガープリント、描画、反応）に整理された、主要な RDKit モジュールと関数の一覧
• references/descriptors_guide.md — 名前、説明、一般的な範囲を含む、200 種類以上の分子記述子の完全な一覧
• references/smarts_patterns.md — 化学クラス別に整理された、官能基検出用の一般的な SMARTS パターン

参考資料

• RDKit Documentation — 公式ドキュメントと入門ガイド
• RDKit Cookbook — レシピと一般的なワークフローパターン
• Getting Started with RDKit in Python — 包括的なチュートリアル
• Lipinski, C.A. et al. (2001) Advanced Drug Delivery Reviews 46:3-26 — ルール・オブ・ファイブ
• Veber, D.F. et al. (2002) J. Med. Chem. 45:2615-2623 — 経口バイオアベイラビリティの基準
• Morgan, H.L. (1965) J. Chem. Doc. 5:107-113 — 環状フィンガープリントアルゴリズム