分子对接模块 · 批量虚拟筛选脚本 vs_batch.py 真实运行记录

验真对象：content/02_分子对接与虚拟筛选/scripts/vs_batch.py（材料自带，声明「仅通过语法检查，未在真实数据上端到端跑过」）
运行环境：Ubuntu 云主机，8 核 / 15 G 内存 / 38 G 可用磁盘，无 docker，micromamba 2.9.0
运行日期：2026-09-23
记录人：小七（子任务执行）
工作目录：verify/docking/（脚本在 run/ 下真实运行；AutoDock-Vina/ 为官方仓库样例来源）

0. 结论速览（先看这段）

问题	结论
脚本是否真跑通	跑通。5 个配体、多进程并行、全部 `ok`，产出 `poses/` 与 `results.csv`
阴性对照	阳性对照伊马替尼重对接得 -13.27 kcal/mol，与官方示例姿态分 -13.234 一致；姿态与共晶结构 RMSD 0.270 Å（官方 0.283 Å），重对接成功
脚本本身是否有问题	有 2 处真实缺陷（同名配体输出互相覆盖 / 输入不存在时不 fail-fast），已修复；原文件存为 `vs_batch.orig.py`，改动见 `diff_vs_batch.md`
跑脚本本身的最小依赖	Python 3 + `vina` 可执行文件 + 输入 PDBQT（脚本只用标准库）；准备输入才需要 meeko / rdkit / molscrub 等
实际耗时	环境创建 ≈ 2 min；单配体试跑 24 s；5 配体批量（exh=8）1 min 17 s；（exh=16）2 min 27 s
如实说明	本记录中所有命令与输出均为真实终端记录，无编造。个别与「预期不同」之处（如 dasatinib 打分低于伊马替尼）在 §6 解释，属打分函数已知局限，非脚本故障

产物清单（均在本目录）： - run_log.md（本文件） - results.csv（脚本产物，final run 输出） - input_manifest.txt（输入文件清单 + md5） - diff_vs_batch.md（脚本改动 diff） - env_versions.txt（关键依赖版本） - pose_rmsd_check.txt（重对接 RMSD 复核） - run/（完整工作目录：input / ligands / poses / 日志 / 边界测试） - AutoDock-Vina/（官方仓库，提供 1IEP 样例数据）

1. 环境搭建

1.1 建独立环境

bash

export MAMBA_ROOT_PREFIX=/root/micromamba_root
micromamba create -y -n dock python=3.11 numpy scipy rdkit vina meeko gemmi openbabel \
    -c conda-forge -c bioconda

结果：Transaction finished，环境建好（日志 create_env.log，耗时约 2 分钟）。

关键版本（完整见 env_versions.txt）：

text

vina 1.2.7           AutoDock Vina f458505-mod
meeko 0.8.0          rdkit 2026.03.1
python 3.11.16       prody 2.6.1
numpy 2.4.6          scipy 1.17.1
gemmi 0.7.5          openbabel 3.2.1

命令行工具就位确认：

bash

$ micromamba run -n dock which vina mk_prepare_receptor.py mk_prepare_ligand.py mk_export.py
/root/micromamba_root/envs/dock/bin/vina
/root/micromamba_root/envs/dock/bin/mk_prepare_receptor.py
/root/micromamba_root/envs/dock/bin/mk_prepare_ligand.py
/root/micromamba_root/envs/dock/bin/mk_export.py

$ micromamba run -n dock vina --version
AutoDock Vina f458505-mod

1.2 按 README 装 molscrub（踩坑见 §5-①）

bash

micromamba run -n dock pip install molscrub      # 报错：缺 joblib
micromamba run -n dock pip install joblib        # 补上后才能跑
micromamba run -n dock scrub.py --help           # OK，显示出 pH / tautomer / ETKDG 等参数

结论：README/requirements.txt 里的 pip install molscrub prody 不完整，scrub.py 会因缺 joblib 直接崩。补 joblib 后正常（建议把 joblib 补进 requirements.txt，见 §7）。

2. 准备真实输入（1IEP 体系）

样例数据来自 AutoDock Vina 官方仓库：

bash

git clone --depth 1 https://github.com/ccsb-scripps/AutoDock-Vina.git
cp AutoDock-Vina/example/basic_docking/data/1iep_receptorH.pdb run/input/
cp AutoDock-Vina/example/basic_docking/data/1iep_ligand.sdf    run/input/

2.1 受体 → PDBQT + 盒子配置

bash

$ micromamba run -n dock mk_prepare_receptor.py -i 1iep_receptorH.pdb -o 1iep_receptor \
      -p -v --box_size 20 20 20 --box_center 15.190 53.903 16.917
@> 4412 atoms and 1 coordinate set(s) were parsed in 0.14s.
Template padding will be used.
 gasteiger charges will be read from template file

Files written:
  1iep_receptor.pdbqt <-- static (i.e., rigid) receptor input file
1iep_receptor.box.txt <-- Vina-style box dimension file
1iep_receptor.box.pdb <-- PDB file to visualize the grid box

生成的盒子配置与官方 solution/1iep_receptor.box.txt 完全一致：

text

center_x = 15.190
center_y = 53.903
center_z = 16.917
size_x = 20.000
size_y = 20.000
size_z = 20.000

2.2 配体：伊马替尼（官方 SDF，共晶坐标）

bash

$ micromamba run -n dock mk_prepare_ligand.py -i 1iep_ligand.sdf -o ../ligands/imatinib.pdbqt
Input molecules processed: 1, skipped: 0
PDBQT files written: 1

产物 3841 字节（与官方 solution/1iep_ligand.pdbqt 同尺寸）。

2.3 配体：再补 4 个（走 SMILES → scrub → SDF → PDBQT 真实流程）

为凑成一张有意义的排序表，除伊马替尼（阳性对照）外补 2 个已知 c-Abl 抑制剂 + 2 个小分子阴性对照：

run/input/extra_ligands.smi：

text

Cc1cn(-c2cc(NC(=O)c3ccc(CN4CCN(C)CC4)cc3C)cc(C(F)(F)F)c2)cn1 nilotinib
Cc1nc(Nc2ncc(C(=O)Nc3c(C)cccc3Cl)s2)cc(N2CCN(CCO)CC2)n1 dasatinib
Cn1c(=O)c2c(ncn2C)n(C)c1=O caffeine
CC(C)Cc1ccc(cc1)C(C)C(=O)O ibuprofen

先用 RDKit 校验 SMILES 合法性，再跑 molscrub 加氢/枚举质子化态/生成 3D：

bash

$ micromamba run -n dock scrub.py extra_ligands.smi -o extra_ligands.sdf --ph 7.4 --cpu 4
Scrub completed.
Input molecules supplied: 4
mols processed: 4, skipped by rdkit: 0, failed: 0
nr isomers (tautomers and acid/base conjugates): 4 (avg. 1.000 per mol)
nr conformers:  4 (avg. 1.000 per isomer, 1.000 per mol)
# real 0m3.543s

从 SDF 批量转 PDBQT：

bash

$ micromamba run -n dock mk_prepare_ligand.py -i extra_ligands.sdf --multimol_outdir ../ligands/
Input molecules processed: 4, skipped: 0
PDBQT files written: 4
No duplicate molecule molecule names were found

最终 run/ligands/ 5 个配体，meeko 在文件头记录了 SMILES，可看到 pH 7.4 下的质子化态：

配体	meeko 写入的 SMILES（节选）	形式电荷	MW	可旋转键
imatinib	`...CN3CC[NH+](C)CC3...`	+1	493.6	7
nilotinib	`...C[N@H+]4CC[N@H+](C)CC4...`	+2	471.5	5
dasatinib	`...N2CC[NH+](CCO)CC2...`	+1	488.0	7
ibuprofen	`...C(=O)[O-]`	−1	206.3	4
caffeine	`Cn1c(=O)c2c(ncn2C)n(C)c1=O`	0	194.2	0

3. 单配体冒烟测试（先确认链路通）

在批量跑之前，先用命令行单独对一次，看打分数量级对不对：

bash

$ time micromamba run -n dock vina --receptor input/1iep_receptor.pdbqt \
      --ligand ligands/imatinib.pdbqt --config input/1iep_receptor.box.txt \
      --exhaustiveness 8 --cpu 4 --seed 42 --out smoke_imatinib_out.pdbqt
Scoring function : vina
Rigid receptor: input/1iep_receptor.pdbqt
Ligand: ligands/imatinib.pdbqt
Grid center: X 15.19 Y 53.903 Z 16.917
Grid size  : X 20 Y 20 Z 20
Grid space : 0.375
Exhaustiveness: 8
CPU: 4
Computing Vina grid ... done.
Performing docking (random seed: 42) ...
mode |   affinity | dist from best mode
     | (kcal/mol) | rmsd l.b.| rmsd u.b.
-----+------------+----------+----------
   1       -13.27          0          0
   2        -11.3      3.045      12.42
...
# real 0m23.607s

结论：-13.27 kcal/mol，落在官方示例「约 −13」的量级。链路通，可以批量跑。

4. 批量对接：用 vs_batch.py 真跑

4.1 原始脚本首次运行（exh=8）

bash

$ python "/Coze/.../scripts/vs_batch.py" \
      --receptor input/1iep_receptor.pdbqt --ligands ligands/ \
      --config input/1iep_receptor.box.txt --out-dir poses --csv results.csv \
      --nproc 5 --cpu-per-job 1 --exhaustiveness 8 --seed 42 --top-n 20

vina        : /root/micromamba_root/envs/dock/bin/vina
受体        : .../run/input/1iep_receptor.pdbqt
config      : .../run/input/1iep_receptor.box.txt
配体数量    : 5
并行进程    : 5（每进程 1 核）
搜索强度    : 8
[1/5] caffeine(ok)
[2/5] ibuprofen(ok)
[3/5] dasatinib(ok)
[4/5] nilotinib(ok)
[5/5] imatinib(ok)

对接结束：成功 5 / 共 5，汇总 -> results.csv

打分最优的前 20 个配体（kcal/mol，越负越可能结合）：
   1. imatinib                         -13.27
   2. nilotinib                         -9.67
   3. dasatinib                         -9.36
   4. ibuprofen                         -8.76
   5. caffeine                          -5.88
# real 1m17.286s

原始脚本在正常输入下是能跑的（这也是它「语法检查通过」之外的真实表现）。

4.2 边界/异常用例探测（找脚本缺陷）

为回答「参数、路径、并发、异常处理」是否有问题，逐个探了 7 个用例：

用例	命令要点	原始脚本表现	判定
T1 `--summarize-only`	只汇总已有 poses	`汇总完成：5 条结果`	正常
T2 `--center/--size` 临时 config	不给 `--config`	生成 `_box.txt` 并跑通 caffeine	正常
T3 不同目录同名配体	`dup1/lig.pdbqt dup2/lig.pdbqt`	报「成功 2 / 共 2」，但 `poses_dup/` 只有 1 个结果、CSV 只 1 行	缺陷，静默丢数据
T4 损坏配体	`garbage` 内容	`broken(failed(rc=1))`，不崩	可接受
T5 受体不存在	`--receptor nope.pdbqt`	逐个配体报 failed，进程 exit 0	缺陷，未 fail-fast
T6 空配体目录	空目录	`没有找到任何配体 PDBQT`，exit 1	正常
T7 vina 路径错	`--vina /no/such/vina`	`caffeine(vina-not-found)`，进程 exit 0	缺陷，未 fail-fast

T3 的关键证据（原始脚本）：

text

[1/2] lig(ok)
[2/2] lig(ok)
对接结束：成功 2 / 共 2，汇总 -> results_dup.csv
--- results_dup.csv ---
rank,ligand,best_affinity_kcal_mol,rmsd_lb,rmsd_ub,n_poses,out_file
1,lig,-8.76,0.0,0.0,7,poses_dup/lig_out.pdbqt      # 只有 1 行，另一个被覆盖

5. 脚本缺陷与修复

原文件已另存为 content/02_分子对接与虚拟筛选/scripts/vs_batch.orig.py（md5 7678bb4f...，与改前 vs_batch.py 完全一致）。完整 diff 见 diff_vs_batch.md，摘要如下。

缺陷 ①：不同目录同名配体，输出互相覆盖（并发/路径）

脚本用「配体文件名去掉扩展名」当输出名，两个都叫 lig.pdbqt 的配体（常见于同时 glob 多个化合物库目录）会写到同一个 poses/lig_out.pdbqt，后跑完的覆盖先跑完的；日志同理会互相冲掉。运行时还按配体数报「成功 N/N」，静默丢数据。

修复：新增 assign_stems()，按 basename 计数；有重名时统一加 __1/__2… 后缀，并在运行前打印提示。单配体或名字本来不重复时，输出命名完全不变（保持 README 里 <配体名>_out.pdbqt 的约定）。run_one() 改为优先用 job 里带上的 stem。

缺陷 ②：输入/依赖缺失时不 fail-fast（异常处理）

受体文件不存在、vina 找不到，都不会在启动前拦住，而是把每个配体都跑一遍、每个都返回 failed，最后 exit 0，容易被误当成"跑完了只是没结果"。

修复：启动前校验 --receptor、--config 文件存在，且 vina 可解析（--vina 路径或 PATH 中可执行），不满足则 argparse.error 直接报清晰中文错误并 exit 2。

修复后复测

bash

# T3 同名配体：不再覆盖，2 个结果都在
提示：以下配体重名，输出名已自动加 __N 后缀避免覆盖：lig
[1/2] lig__1(ok)
[2/2] lig__2(ok)
对接结束：成功 2 / 共 2，汇总 -> results_dup.csv
--- results_dup.csv ---
rank,ligand,best_affinity_kcal_mol,rmsd_lb,rmsd_ub,n_poses,out_file
1,lig__2,-8.76,0.0,0.0,7,poses_dup/lig__2_out.pdbqt
2,lig__1,-5.878,0.0,0.0,9,poses_dup/lig__1_out.pdbqt

# T5 受体不存在：直接拦下
vs_batch.py: error: 找不到受体文件：input/nope.pdbqt        (exit 2)

# T7 vina 路径错：直接拦下
vs_batch.py: error: 找不到 vina 可执行文件（/no/such/vina）。请先 conda install -c conda-forge vina，或用 --vina 指定完整路径。   (exit 2)

其余用例（T1/T2/T4/T6）修复后行为不变。语法与 import 均通过 python -m py_compile。

6. final run 与结果合理性

用修复后的脚本做一次正式批量（试跑档 exh=16，固定种子）：

bash

$ time python ".../vs_batch.py" --receptor input/1iep_receptor.pdbqt --ligands ligands/ \
      --config input/1iep_receptor.box.txt --out-dir poses --csv results.csv \
      --nproc 5 --cpu-per-job 1 --exhaustiveness 16 --seed 42 --top-n 20

配体数量    : 5
并行进程    : 5（每进程 1 核）
搜索强度    : 16
[1/5] caffeine(ok)   [2/5] ibuprofen(ok)   [3/5] dasatinib(ok)
[4/5] nilotinib(ok)  [5/5] imatinib(ok)
对接结束：成功 5 / 共 5，汇总 -> results.csv
   1. imatinib                         -13.27
   2. nilotinib                         -9.67
   3. dasatinib                         -9.36
   4. ibuprofen                         -8.76
   5. caffeine                          -5.86
# real 2m27.089s

results.csv（= 本目录 results.csv）：

text

rank,ligand,best_affinity_kcal_mol,rmsd_lb,rmsd_ub,n_poses,out_file
1,imatinib,-13.273,0.0,0.0,4,poses/imatinib_out.pdbqt
2,nilotinib,-9.666,0.0,0.0,9,poses/nilotinib_out.pdbqt
3,dasatinib,-9.363,0.0,0.0,9,poses/dasatinib_out.pdbqt
4,ibuprofen,-8.76,0.0,0.0,9,poses/ibuprofen_out.pdbqt
5,caffeine,-5.86,0.0,0.0,9,poses/caffeine_out.pdbqt

6.1 阳性对照复核（打分 + 姿态双重）

我们自己跑的伊马替尼 top 姿态：-13.273
官方示例 solution/1iep_ligand_vina_out.pdbqt top：-13.234
两者姿态与 1IEP 共晶配体的重原子 RMSD（rdMolAlign.GetBestRMS，对称性校正，见 pose_rmsd_check.txt）：

text

crystal ligand: input/1iep_ligand.sdf (1IEP 共晶伊马替尼), 重原子数 = 37
our poses: 4  poses;  official example poses: 4  poses
RMSD(our imatinib top pose  vs crystal) = 0.270 A
RMSD(official top pose       vs crystal) = 0.283 A

打分与官方一致、姿态 RMSD < 0.3 Å（远好于常用的 < 2 Å 判据）——重对接成功，说明从受体准备到批量脚本整条链路是对的。

6.2 关于排序里「反直觉」的一点（如实说明，非 bug）

nilotinib / dasatinib 是公认的 c-Abl 抑制剂，分数（-9.7 / -9.4）却低于伊马替尼。原因不在脚本：

Vina 打分函数对难分伯仲的同类抑制剂排序能力本就有限，绝对值不能当 Kd（README 第四节自己也写了）；
更大的影响来自 pH 7.4 质子化态——scrub 把 nilotinib 变成双电荷（+2）、dasatinib 单电荷（+1），Vina 对带电配体在疏水口袋里的去溶剂化惩罚会明显压低分数；伊马替尼只有 +1 且更贴合共晶姿态；
盒子里伊马替尼本来就是共晶分子，geometry 天生占优。

想让排序更靠谱，按 README/笔记的建议：叠加关键残基相互作用过滤、挑候选做 MM-GBSA 或短 MD 复算。这不是脚本要修的东西。

7. 给主 agent / 内容团队的结论与建议

脚本判定：真实可跑通，可上网站。带 5 配体（含阳性对照）批量一次成功，产物齐全。
必改的脚本缺陷已改（同名覆盖、缺 fail-fast），原文件保留为 vs_batch.orig.py，diff 见 diff_vs_batch.md。
文档还有 1 处小坑（不在授权改动范围内，仅建议）：requirements.txt 与 README.md 的 pip install molscrub prody 漏了 joblib，照此安装后 scrub.py 会 ModuleNotFoundError: No module named 'joblib'。建议在 requirements.txt 的 molscrub 一行附近补 joblib，README 的 pip 命令改成 pip install molscrub joblib prody。
最小依赖： - 只跑 vs_batch.py：Python 3 + vina 可执行（脚本仅用标准库）+ 准备好的 PDBQT/config。 - 从原始结构准备输入：受体 mk_prepare_receptor.py（meeko，依赖 prody）；配体 mk_prepare_ligand.py（meeko）或 SMILES 路线 scrub.py（molscrub + joblib + rdkit）。
给初学者的操作要点： - 先单独跑 1 个配体冒烟测试，确认打分在合理量级，再上批量； - 一定带阳性对照（本项目里就是伊马替尼），否则结果无从判断； - --out-dir 别和输入目录混放；--nproc × --cpu-per-job 别超过物理核数（本机 8 核，5×1 合适）； - 固定 --seed 保证可复现；中断后用 --summarize-only 重新汇总； - 配体重名会被加 __N 后缀（修复后），看到后缀属正常； - --exhaustiveness 试跑用 8–16，正式筛选 32，别一上来就拉满。

附：本记录对应的关键命令一览

bash

export MAMBA_ROOT_PREFIX=/root/micromamba_root
micromamba create -y -n dock python=3.11 numpy scipy rdkit vina meeko gemmi openbabel -c conda-forge -c bioconda
micromamba run -n dock pip install molscrub joblib
git clone --depth 1 https://github.com/ccsb-scripps/AutoDock-Vina.git
micromamba run -n dock mk_prepare_receptor.py -i 1iep_receptorH.pdb -o 1iep_receptor -p -v --box_size 20 20 20 --box_center 15.190 53.903 16.917
micromamba run -n dock mk_prepare_ligand.py -i 1iep_ligand.sdf -o ../ligands/imatinib.pdbqt
micromamba run -n dock scrub.py extra_ligands.smi -o extra_ligands.sdf --ph 7.4 --cpu 4
micromamba run -n dock mk_prepare_ligand.py -i extra_ligands.sdf --multimol_outdir ../ligands/
python scripts/vs_batch.py --receptor input/1iep_receptor.pdbqt --ligands ligands/ --config input/1iep_receptor.box.txt --out-dir poses --csv results.csv --nproc 5 --cpu-per-job 1 --exhaustiveness 16 --seed 42 --top-n 20

（所有输出均来自以上命令的真实执行，未作删改。）