这篇把对接流程串成一条可直接运行的链路：结构准备 → 重对接验证 → 批量筛选 → 结果分析。其中重对接验证是绝不能跳过的一步——它花几分钟，却能挡掉大部分「跑完了但全是错的」情况。
环境准备
pip install vina meeko rdkit
mamba install -c conda-forge pdbfixer openmm # 结构修复
# 或直接用 GNINA 二进制（能读 PDB/SDF，省掉 PDBQT 转换）第一步：受体准备
from pdbfixer import PDBFixer
from openmm.app import PDBFile
fixer = PDBFixer(pdbid="1M17") # EGFR 激酶域 + 厄洛替尼类似物
fixer.findMissingResidues()
print("缺失残基:", fixer.missingResidues)
fixer.findNonstandardResidues()
fixer.replaceNonstandardResidues()
fixer.removeHeterogens(keepWater=False) # 注意：这会删掉配体
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(pH=7.4)
PDBFile.writeFile(fixer.topology, fixer.positions,
open("receptor_fixed.pdb", "w"), keepIds=True)# 转成 PDBQT
mk_prepare_receptor.py -i receptor_fixed.pdb -o receptor -p注意 removeHeterogens 会删掉共晶配体——所以要先把配体单独提取出来，用于定盒子和重对接验证。
第二步：确定对接盒子
from rdkit import Chem
import numpy as np
# 从原始 PDB 中提取共晶配体（用 PyMOL 或手工分离）
ref = Chem.MolFromMolFile("ref_ligand.sdf", removeHs=False)
conf = ref.GetConformer()
coords = np.array([list(conf.GetAtomPosition(i))
for i in range(ref.GetNumAtoms())])
center = coords.mean(axis=0)
extent = coords.max(axis=0) - coords.min(axis=0)
box = extent + 10.0 # 配体最长径 + 10 Å 余量
box = np.maximum(box, 20.0) # 最小 20 Å
print(f"中心: {center.round(2)}")
print(f"盒子: {box.round(1)}")盒子设置是最高频的错误来源：太小切掉正确姿势，太大稀释搜索、增加噪声。有共晶配体时直接用它的质心最可靠。
第三步：重对接验证（必做）
from vina import Vina
from rdkit.Chem import AllChem, rdMolAlign
# 准备参考配体的 PDBQT
# mk_prepare_ligand.py -i ref_ligand.sdf -o ref_ligand.pdbqt
v = Vina(sf_name="vina", seed=42, cpu=8)
v.set_receptor("receptor.pdbqt")
v.compute_vina_maps(center=center.tolist(), box_size=box.tolist())
v.set_ligand_from_file("ref_ligand.pdbqt")
v.dock(exhaustiveness=32, n_poses=10)
v.write_poses("redock_out.pdbqt", n_poses=10, overwrite=True)
# 转回 SDF 并算 RMSD
# mk_export.py redock_out.pdbqt -s redock.sdf
docked = Chem.SDMolSupplier("redock.sdf", removeHs=True)[0]
ref_noh = Chem.RemoveHs(ref)
rmsd = rdMolAlign.CalcRMS(Chem.RemoveHs(docked), ref_noh)
print(f"重对接 RMSD = {rmsd:.2f} Å")
if rmsd <= 2.0:
print("✓ 流程验证通过，可以做新分子")
else:
print("✗ 未复现晶体姿势 —— 检查结构准备、质子化态、盒子设置")
print(" 在解决之前，对新分子的结果没有意义")判据是重原子 RMSD ≤ 2.0 Å。做不到就说明准备或盒子有问题，此时对新分子的对接结果完全不可信。
第四步：批量筛选
import glob, json
from vina import Vina
from concurrent.futures import ProcessPoolExecutor
def dock_one(args):
pdbqt_path, center, box = args
try:
v = Vina(sf_name="vina", seed=42, cpu=1, verbosity=0)
v.set_receptor("receptor.pdbqt")
v.compute_vina_maps(center=center, box_size=box)
v.set_ligand_from_file(pdbqt_path)
v.dock(exhaustiveness=8, n_poses=5) # 初筛用 8 即可
energies = v.energies()
name = pdbqt_path.split("/")[-1].replace(".pdbqt", "")
return {"name": name, "score": float(energies[0][0])}
except Exception as e:
return {"name": pdbqt_path, "score": None, "error": str(e)}
ligands = glob.glob("ligands_pdbqt/*.pdbqt")
args = [(p, center.tolist(), box.tolist()) for p in ligands]
with ProcessPoolExecutor(max_workers=16) as ex:
results = list(ex.map(dock_one, args, chunksize=10))
import pandas as pd
df = pd.DataFrame([r for r in results if r.get("score") is not None])
df = df.sort_values("score")
print(f"成功 {len(df)}/{len(ligands)}；最佳分数 {df['score'].min():.2f}")
df.to_csv("docking_results.csv", index=False)第五步：结果分析
# 1) 配体高效性 —— 比原始分数公平
from rdkit.Chem import Descriptors
df["n_heavy"] = df["smiles"].map(
lambda s: Chem.MolFromSmiles(s).GetNumHeavyAtoms())
df["LE"] = -df["score"] / df["n_heavy"] # 每重原子的结合贡献
# Vina 分数天然偏好大分子；LE 能纠正这个偏倚
print(df.nlargest(10, "LE")[["name", "score", "n_heavy", "LE"]])
# 2) 分数分布 —— 检查是否合理
print(df["score"].describe())
# 常见活性分子在 -7 ~ -11；若大量分子低于 -13，可能盒子有问题
# 3) 与已知活性分子对比（如果有）
# 已知活性分子应该排在前面，否则流程有问题怎么读 Vina 分数
- 单位是 kcal/mol 的估算结合自由能，越负越好；
- 绝对值不可当亲和力用——与实测 IC50 的相关性通常很弱（Pearson 常在 0.3～0.5）；
- 评分函数偏好大分子和高疏水性分子，比较不同大小的分子时看配体高效性更公平；
- 正确用法是粗筛排序 + 姿势假设生成，不是排名预测。
常见坑与提示
- 重对接验证（RMSD ≤ 2 Å）是必做的第一步，不通过就别做新分子；
- 盒子中心用共晶配体质心，尺寸 = 配体最长径 + 10 Å；
- 固定
seed与exhaustiveness，报告时一并写明； - 比较不同大小分子时用配体高效性，别直接比原始分数。