diffdock

diffdock

概述

DiffDock 使用扩散生成模型，直接根据蛋白质结构和配体 SMILES 预测蛋白质-配体结合构象，将对接视为生成问题而非搜索问题。与传统对接工具（AutoDock Vina、Glide）不同，DiffDock 不需要预先指定结合位点，而是在整个蛋白质表面采样构象。它会输出一组经过排序的结合构象及其对应的置信度分数。DiffDock 擅长盲对接任务，并能生成多样的构象假设，因此适用于从头发现结合位点以及处理具有挑战性的靶标。

适用场景

• 盲对接（未知结合位点）：不知道配体会与蛋白质的哪个位置结合，希望发现候选结合位点。
• 传统对接难以处理的靶标：变构位点、柔性区域，或靶标结合位点没有共晶结构的蛋白质。
• 探索多种结合模式：生成多样的构象集合，以了解结合事件中的构象灵活性。
• 构效关系（SAR）探索：快速对接一系列类似物，以比较预测的结合模式。
• 生成片段筛选假设：为片段分子识别合理的结合位点。
• 对于结合位点已知且假设蛋白质为刚性的情况，AutoDock Vina 或 GNINA 可能速度更快，准确度也相当。
• 对于大规模虚拟筛选（>10,000 个化合物），建议使用 GNINA 或 DiffDock-L（大规模版本），而非标准 DiffDock。
• 当结合口袋定义明确，且大型化合物库需要更高处理速度时，请改用 AutoDock Vina

前提条件

• Python 包：diffdock（建议通过 conda 安装）、rdkit、torch、biopython、nglview（可视化）
• 系统：强烈建议使用 GPU（NVIDIA CUDA）；CPU 推理速度较慢（每个化合物约需 5–10 分钟）
• 数据要求：蛋白质 PDB 文件（已清理并添加质子），配体为 SMILES 字符串或 SDF 文件
• 环境：包含兼容 CUDA 的 PyTorch 的 conda 环境

# Recommended: clone and install from source
git clone https://github.com/gcorso/DiffDock.git
cd DiffDock
conda create -n diffdock python=3.9
conda activate diffdock
pip install torch torchvision --extra-index-url https://download.pytorch.org/whl/cu118
pip install -r requirements.txt

# Download pretrained model weights
python -c "from utils.download import download_pretrained; download_pretrained()"

工作流

步骤 1：准备蛋白质结构

from Bio import PDB
from Bio.PDB import PDBParser, PDBIO, Select

class NonHetSelect(Select):
    """Remove HETATM records (ligands, water) — keep only protein atoms."""
    def accept_residue(self, residue):
        return residue.id[0] == " "

def clean_pdb(input_pdb: str, output_pdb: str):
    parser = PDBParser(QUIET=True)
    structure = parser.get_structure("protein", input_pdb)
    io = PDBIO()
    io.set_structure(structure)
    io.save(output_pdb, NonHetSelect())
    print(f"Cleaned PDB saved to: {output_pdb}")

clean_pdb("raw_protein.pdb", "protein_clean.pdb")

步骤 2：准备配体输入

from rdkit import Chem
from rdkit.Chem import AllChem, SDWriter

def smiles_to_sdf(smiles: str, output_sdf: str, n_confs: int = 1):
    """Convert SMILES to 3D SDF for DiffDock input."""
    mol = Chem.MolFromSmiles(smiles)
    mol = Chem.AddHs(mol)
    AllChem.EmbedMolecule(mol, AllChem.ETKDGv3())
    AllChem.MMFFOptimizeMolecule(mol)
    writer = SDWriter(output_sdf)
    writer.write(mol)
    writer.close()
    print(f"Ligand SDF written to: {output_sdf}")
    return mol

# Example: ibuprofen
smiles = "CC(C)Cc1ccc(cc1)C(C)C(=O)O"
mol    = smiles_to_sdf(smiles, "ligand.sdf")
print(f"Ligand formula: {Chem.rdMolDescriptors.CalcMolFormula(mol)}")

步骤 3：运行 DiffDock 推理

# Command-line inference (run from the DiffDock directory)
python inference.py \
    --protein_path protein_clean.pdb \
    --ligand       "CC(C)Cc1ccc(cc1)C(C)C(=O)O" \
    --out_dir      results/ \
    --inference_steps 20 \
    --samples_per_complex 40 \
    --batch_size 10 \
    --no_final_step_noise

import subprocess

def run_diffdock(protein_pdb: str, ligand_smiles: str, out_dir: str,
                 n_samples: int = 40, n_steps: int = 20):
    cmd = [
        "python", "inference.py",
        "--protein_path",      protein_pdb,
        "--ligand",            ligand_smiles,
        "--out_dir",           out_dir,
        "--inference_steps",   str(n_steps),
        "--samples_per_complex", str(n_samples),
        "--batch_size",        "10",
        "--no_final_step_noise",
    ]
    result = subprocess.run(cmd, capture_output=True, text=True, cwd="DiffDock/")
    if result.returncode == 0:
        print(f"DiffDock complete. Results in: {out_dir}")
    else:
        print(f"Error: {result.stderr}")
    return result

run_diffdock("protein_clean.pdb", "CC(C)Cc1ccc(cc1)C(C)C(=O)O", "results/")

步骤 4：解析并排列置信度分数

import re
from pathlib import Path
import pandas as pd

def parse_diffdock_results(out_dir: str) -> pd.DataFrame:
    """Parse DiffDock output SDF files and confidence scores."""
    out_path = Path(out_dir)
    records  = []

    # DiffDock names output files: rank{N}_confidence{score}.sdf
    for sdf_file in sorted(out_path.glob("rank*_confidence*.sdf")):
        name = sdf_file.stem
        # Extract rank and confidence from filename
        rank_match = re.search(r"rank(\d+)", name)
        conf_match = re.search(r"confidence(-?[\d.]+)", name)
        if rank_match and conf_match:
            records.append({
                "rank":       int(rank_match.group(1)),
                "confidence": float(conf_match.group(1)),
                "sdf_file":   str(sdf_file),
            })

    df = pd.DataFrame(records).sort_values("rank")
    print(f"Found {len(df)} poses")
    print(df[["rank", "confidence", "sdf_file"]].head(10))
    return df

df_results = parse_diffdock_results("results/")

步骤 5：分析排名靠前的构象 — 提取结合位点残基

from rdkit import Chem
from rdkit.Chem import AllChem
from Bio.PDB import PDBParser
import numpy as np

def get_binding_residues(protein_pdb: str, ligand_sdf: str, cutoff_angstrom: float = 4.0):
    """Find protein residues within cutoff distance of the top-ranked ligand pose."""
    parser    = PDBParser(QUIET=True)
    structure = parser.get_structure("prot", protein_pdb)
    prot_atoms = [(atom.get_coord(), residue.resname, residue.id[1])
                  for chain in structure for residue in chain
                  for atom in residue.get_atoms()]

    mol = Chem.SDMolSupplier(ligand_sdf, removeHs=False)[0]
    lig_coords = mol.GetConformer().GetPositions()

    contacts = []
    for prot_coord, resname, resnum in prot_atoms:
        dists = np.linalg.norm(lig_coords - prot_coord, axis=1)
        if dists.min() <= cutoff_angstrom:
            contacts.append((resnum, resname))

    contacts = sorted(set(contacts))
    print(f"Binding site residues within {cutoff_angstrom} A: {contacts[:10]}")
    return contacts

# Use top-ranked pose
top_sdf = df_results.loc[df_results.rank == 1, "sdf_file"].iloc[0]
contacts = get_binding_residues("protein_clean.pdb", top_sdf)

步骤 6：在 NGLview 中可视化构象

import nglview as nv
from rdkit import Chem

# Load protein + top pose in Jupyter notebook
view = nv.NGLWidget()
view.add_pdbfile("protein_clean.pdb")

top_sdf  = df_results.loc[df_results.rank == 1, "sdf_file"].iloc[0]
view.add_component(top_sdf)
view.representations = [
    {"type": "cartoon", "params": {"color": "chainindex"}},
    {"type": "ball+stick", "params": {"sele": "ligand"}},
]
print(f"Visualizing top pose: confidence={df_results.confidence.iloc[0]:.3f}")
view

关键参数

参数	默认值	范围 / 选项	作用
--inference_steps	20	10–40	扩散反向步骤数；步数越多，速度越慢，但准确度越高
--samples_per_complex	40	10–100	采样的构象数；数量越多，对结合模式的覆盖越好
--batch_size	10	1–32	GPU 批次大小；发生 OOM 错误时可调低
--no_final_step_noise	off	flag	移除扩散最后一步的噪声；提升构象质量
--actual_steps	equals inference_steps	1–inference_steps	实际运行的步数（可以少于总步数）
--save_visualisation	off	flag	除 SDF 外，还保存 PDB 可视化文件
cutoff_angstrom	4.0	3.0–6.0 Å	用于定义结合位点残基的距离截断值

常用示例

示例：批量对接多个配体

适用场景：将一组类似物对接到同一种蛋白质，以进行 SAR 分析。

import pandas as pd
import subprocess

smiles_list = [
    ("compound_1", "CC(C)Cc1ccc(cc1)C(C)C(=O)O"),
    ("compound_2", "CC(C)Cc1ccc(cc1)C(C)C(=O)N"),
    ("compound_3", "CC(C)Cc1ccc(cc1)C(C)C(=O)OC"),
]

results = []
for name, smiles in smiles_list:
    out = f"results/{name}"
    cmd = ["python", "inference.py",
           "--protein_path", "protein_clean.pdb",
           "--ligand", smiles,
           "--out_dir", out,
           "--inference_steps", "20",
           "--samples_per_complex", "20"]
    subprocess.run(cmd, cwd="DiffDock/", capture_output=True)
    # Parse top confidence score
    df_r = parse_diffdock_results(out)
    if not df_r.empty:
        top_conf = df_r.loc[df_r.rank == 1, "confidence"].iloc[0]
        results.append({"name": name, "smiles": smiles, "top_confidence": top_conf})

df_batch = pd.DataFrame(results).sort_values("top_confidence", ascending=False)
df_batch.to_csv("batch_docking_results.csv", index=False)
print(df_batch)

示例：按置信度阈值筛选构象

适用场景：仅保留高置信度构象，以便进一步分析或可视化。

# Confidence > 0 generally indicates a plausible binding pose
# DiffDock confidence scores: higher = more confident; ~0 is marginal; < -1 is poor
high_conf = df_results[df_results["confidence"] > 0.0]
print(f"High-confidence poses: {len(high_conf)} / {len(df_results)}")
print(high_conf[["rank", "confidence", "sdf_file"]])

示例：将排名最高的构象转换为 PDBQT，以便使用 Vina 重新评分

适用场景：使用 AutoDock Vina 的能量函数对 DiffDock 构象重新评分。

# Convert SDF to PDBQT using OpenBabel
obabel rank1_confidence0.75.sdf -O rank1_ligand.pdbqt
obabel protein_clean.pdb -O protein.pdbqt -xr

# Rescore (no docking search, just energy evaluation)
vina --receptor protein.pdbqt --ligand rank1_ligand.pdbqt \
     --score_only --out rank1_rescored.pdbqt

预期输出

• results/rank{N}_confidence{score}.sdf — 按置信度得分排序的 3D 配体构象
• df_results DataFrame，包含每个构象的排名、置信度得分和文件路径
• 结合位点残基列表（截止距离范围内的残基编号和名称）
• Jupyter 笔记本中的 NGLview 交互式 3D 可视化

故障排除

问题	原因	解决方案
CUDA out of memory	批次大小对于 GPU 而言过大	将 --batch_size 减少到 4 或 2
结果目录为空	蛋白质 PDB 解析失败	确保 PDB 仅包含 ATOM 记录；使用 clean_pdb() 移除 HETATM
所有置信度得分均 < -2	配体或蛋白质格式问题	使用 RDKit 验证 SMILES；确保蛋白质已质子化且完整
推理速度极慢（>30 min）	正在 CPU 上运行	强烈建议使用 GPU；必须正确配置 CUDA 环境
ModuleNotFoundError: e3nn	未安装依赖项	在 DiffDock conda 环境中运行 pip install e3nn
构象聚集在同一位点	--samples_per_complex 值过低	增大到 40–100，以覆盖更多位点
蛋白质缺少残基	晶体结构不完整	在对接前使用 MODELLER 或 Swiss-Model 补齐缺失部分

参考资料

• DiffDock GitHub（gcorso/DiffDock） — 源代码、安装说明和预训练权重
• Corso 等人（2023），ICLR — DiffDock 论文 — 原始扩散对接论文
• DiffDock Web 服务 — 无需本地安装即可通过浏览器进行推理
• PatentsView API 文档 — 用于知识产权分析