文章/答案/技术大牛

发布

社区首页 >专栏 >课前准备--采用pyrosetta对autodock vina的分子对接结果进行优化

课前准备--采用pyrosetta对autodock vina的分子对接结果进行优化

原创

追风少年i

发布于 2026-08-15 11:12:18

1280

作者，Evil Genius

接下来要为我们的生化课程，分子对接/分子动力学做好准备，大家学习分子对接/分子动力学，autodock vina、pyrosetta、Gromacs是必要要掌握的。

这一篇我们来采用pyrosetta对autodock vina的分子对接结果进行优化/

首先了解一下背景

Vina 不是简单的“刚性对接”。

更准确地说，经典 AutoDock Vina 是 rigid receptor + flexible ligand：受体通常保持固定，而 ligand 的平移、旋转和可旋转键构象可以搜索。因此，Vina 最大的局限之一不是“配体不能动”，而是蛋白结合口袋的柔性没有被充分考虑。Vina 官方文档也将 receptor rigidity 与 ligand flexibility 分开处理。

一、首先理解：Vina 到底在做什么？

假设蛋白：

Protein
│
└── Binding pocket
       │
       ├── Phe
       ├── Tyr
       ├── Asp
       ├── Leu
       └── Arg

ligand：

LIG

Vina 要解决的问题是：

在指定 binding box 中，寻找 ligand 比较合理的空间位置和构象。

可以理解为：

Ligand
  │
  ├── translation
  ├── rotation
  └── torsion
       │
       ▼
Binding pocket
       │
       ▼
不同 docking poses

最后得到：

Pose 1    -10.2 kcal/mol
Pose 2     -9.8 kcal/mol
Pose 3     -9.5 kcal/mol
...

二、Vina 的第一个核心局限：受体通常是刚性的

假设真实结合过程：

未结合：

       Phe
        \
         \
          pocket

Ligand进入：

真实情况下：

ligand 进入 binding pocket 后，蛋白侧链可能发生 rotamer rearrangement。

但是经典 Vina：

Phe ── 固定
Tyr ── 固定
Asp ── 固定

        ↓

      LIG
        ↓
寻找最佳位置

也就是说：

Protein
├── backbone       固定
├── side chain     基本固定
│
└── ligand
      ├── rotation
      ├── translation
      └── torsion

所以它更准确的定义是：

rigid-receptor / flexible-ligand docking

而不是：

rigid ligand docking。

三、第二个局限：Vina 的 binding pocket 不会充分“诱导适配”

蛋白-配体结合经常存在：

induced fit

例如：

Vina：

PHE100
   \
    \
     LIG


真实：

PHE100
   ↻
    \
     LIG

也就是：

Ligand
   ↓
接近 binding pocket
   ↓
side-chain rearrangement
   ↓
更合理的 interface

Vina 的普通 rigid receptor 模式不能充分描述这一过程。

因此：

Vina 得到的是一个 docking pose，而不一定是最终能量最低的 protein-ligand complex。

四、第三个局限：Vina 的最佳 pose 不一定是真实 pose

例如：

Pose 1   -10.2
Pose 2    -9.9
Pose 3    -9.7

不能直接认为：

Pose 1 = 正确
Pose 2 = 错误

因为 docking score 是一个近似评分函数。

它并不是严格的：

实验结合自由能

而是：

近似 scoring function

所以可能出现：

Pose A

Vina = -10.5

但是：

steric clash
+
不合理 side chain
+
H-bond geometry不好

而：

Pose B

Vina = -9.8

但是：

H-bond合理
+
interface packing好
+
没有明显clash

因此：

不能只看 Vina score。

五、第四个局限：蛋白侧链构象可能来自 apo 状态

比如蛋白 PDB 是：

apo protein

那么：

Phe100
Tyr105
Asp110

的构象可能是：

apo state

但是 ligand 结合后：

holo state

可能不同。

于是：

apo protein
      ↓
Vina
      ↓
ligand docking

得到的结果可能受到：

初始蛋白构象的限制。

六、第五个局限：Vina 不是 MD

Vina 做的是：

搜索 + 打分

而不是：

时间尺度上的动力学模拟

所以它不会告诉：

这个 ligand 在 100 ns 后是否稳定？

binding pocket 是否持续存在？

H-bond 是否持续？

protein backbone 是否发生变化？

这些问题应该由：

GROMACS / AMBER 等 MD

进一步解决。

七、那么为什么使用 PyRosetta？

这里就出现了一个非常自然的逻辑：

Vina
│
│  ligand pose search
▼
初始 complex
│
│  但 protein side chain 没有充分优化
▼
PyRosetta
│
├── side-chain repacking
├── local minimization
├── protein-ligand interface optimization
└── local relaxation
│
▼
refined complex
│
▼
GROMACS
│
└── dynamic stability

所以：

Vina 负责“找到一个合理的结合姿势”，PyRosetta 负责“在这个姿势附近把蛋白-配体界面调整得更合理”。

八、为什么一定强调“局部优化”？

因为并不是想重新折叠整个蛋白。

假设：

Protein = 500 aa

ligand 只和：

20 residues

发生直接相互作用。

那么真正关心的是：

Ligand
  ↓
5~8 Å
  ↓
Binding pocket

而不是：

整个 500 aa protein

因此：

远离 ligand 的区域
        ↓
尽量固定

ligand 周围 5~8 Å
        ↓
允许调整

这就是：

local refinement

九、完整工作流程

最终推荐：

                    Protein PDB
                         │
                         ▼
                 Protein preparation
                         │
                         ▼
                    AutoDock Vina
                         │
              ┌──────────┼──────────┐
              ▼          ▼          ▼
            Pose1      Pose2      Pose3
              │          │          │
              └──────────┼──────────┘
                         ▼
                  Top 5~10 poses
                         │
                         ▼
                Ligand parametrization
                         │
                         ▼
                    LIG.params
                         │
                         ▼
              Protein + ligand complex
                         │
                         ▼
                  PyRosetta Pose
                         │
                         ▼
              Identify 5~8 Å pocket
                         │
                         ▼
                  Side-chain repack
                         │
                         ▼
                  Local minimization
                         │
                         ▼
                    Local Relax
                         │
                         ▼
                  Final minimization
                         │
                         ▼
                 Refined complexes
                         │
                         ▼
             Interface / H-bond / clash
                         │
                         ▼
                   Pose ranking
                         │
                         ▼
                      GROMACS

十、Step 1：准备文件

建议目录：

vina_rosetta/
├── input/
│   ├── protein.pdb
│   └── ligand.sdf
│
├── vina/
│   ├── receptor.pdbqt
│   ├── ligand.pdbqt
│   ├── config.txt
│   └── vina_out.pdbqt
│
├── params/
│   ├── LIG.params
│   └── LIG_0001.pdb
│
├── poses/
│
├── complexes/
│
├── optimized/
│
└── results/

十一、Step 2：Vina docking

例如：

vina \
    --receptor vina/receptor.pdbqt \
    --ligand vina/ligand.pdbqt \
    --config vina/config.txt \
    --exhaustiveness 32 \
    --num_modes 10 \
    --out vina/vina_out.pdbqt

这里：

--num_modes 10

让 Vina 输出多个 pose。

建议：

Top 5

作为最低限度。

更稳妥：

Top 10

十二、Step 3：为什么保留 Top 10？

因为：

Vina score

不是绝对真值。

例如：

Pose 1   -10.4
Pose 2   -10.1
Pose 3   -9.9
Pose 4   -9.8
Pose 5   -9.6

这些 pose 可能代表：

同一个 binding mode 的不同构象

也可能代表：

不同 binding modes

所以后续 PyRosetta 应该：

Pose 1 → refine
Pose 2 → refine
Pose 3 → refine
...

而不是：

只 refine Pose 1

十三、Step 4：准备 ligand 参数

这是 PyRosetta 最重要的准备工作之一。

Vina 使用：

ligand.pdbqt

PyRosetta 使用：

LIG.params

所以：

ligand.sdf
     │
     ▼
molfile_to_params.py
     │
     ├── LIG.params
     └── LIG_0001.pdb

Rosetta 官方的 ligand preparation 流程就是通过 molfile_to_params.py 从分子文件生成 ligand params。

例如：

python \
$ROSETTA/main/source/scripts/python/public/molfile_to_params.py \
-n LIG \
-p LIG \
input/ligand.sdf

得到：

LIG.params
LIG_0001.pdb

十四、.params 到底是什么？

可以把它理解成：

告诉 Rosetta：“LIG 是什么东西”。

包括：

原子
│
├── C
├── N
├── O
└── ...

键
│
├── C-C
├── C-N
└── C-O

torsion
│
└── 哪些键可以旋转

Rosetta atom type

因此：

LIG.pdb

只有坐标是不够的。

Rosetta 还需要：

LIG.params

来理解 ligand 的化学结构。

十五、Step 5：Vina pose 和 Rosetta ligand 坐标怎么结合？

这里是整个 pipeline 中最容易出错的部分。

实际上有两个信息来源：

化学拓扑

来自：

ligand.sdf

↓

LIG.params

空间坐标

来自：

Vina pose

所以：

SDF
 ↓
化学结构

Vina
 ↓
空间位置

最终：

Protein
 +
Vina ligand coordinates
 +
Rosetta ligand topology

组成：

complex.pdb

不要简单地把 PDBQT 改成 PDB 后就结束。

十六、Step 6：PyRosetta 初始化

import pyrosetta

pyrosetta.init(
    "-extra_res_fa params/LIG.params "
    "-mute all"
)

然后：

from pyrosetta import *
from pyrosetta.rosetta import *

加载：

pose = pose_from_pdb(
    "complexes/complex_01.pdb"
)

检查：

for i in range(
    1,
    pose.total_residue() + 1
):
    print(
        i,
        pose.residue(i).name3()
    )

应该看到：

1 ALA
2 GLY
3 VAL
...
250 TYR
251 LIG

十七、Step 7：找到 ligand

ligand_res = None

for i in range(
    1,
    pose.total_residue() + 1
):

    if pose.residue(i).name3() == "LIG":

        ligand_res = i
        break

if ligand_res is None:
    raise RuntimeError(
        "LIG not found"
    )

print(
    "Ligand residue:",
    ligand_res
)

十八、Step 8：找 binding pocket

假设使用：

6 Å

作为局部优化范围。

from pyrosetta.rosetta.core.select import residue_selector

lig_selector = (
    residue_selector.ResidueIndexSelector(
        str(ligand_res)
    )
)

pocket_selector = (
    residue_selector.NeighborhoodResidueSelector()
)

pocket_selector.set_focus(
    lig_selector
)

pocket_selector.set_distance(
    6.0
)

pocket_selector.set_include_focus_in_subset(
    True
)

得到 pocket：

subset = pocket_selector.apply(
    pose
)

pocket_residues = []

for i in range(
    1,
    pose.total_residue() + 1
):

    if subset[i]:

        pocket_residues.append(i)

print(
    pocket_residues
)

十九、为什么是 6 Å？

不是一个绝对规则。

可以尝试：

4 Å
5 Å
6 Å
8 Å

通常：

4 Å

更严格。

6 Å

是一个比较实用的起点。

8 Å

允许更大范围的局部适配。

对于第一次建立 pipeline：

二十、Step 9：计算初始 Rosetta score

scorefxn = get_fa_scorefxn()

initial_score = scorefxn(
    pose
)

print(
    "Initial score:",
    initial_score
)

注意：

这个 score 和 Vina score 不能直接比较。

例如：

Vina = -10.2
Rosetta = -135

不能说：

-135 比 -10.2 好

它们是完全不同的 scoring function。

二十一、Step 10：Side-chain Repacking

这是第一个真正的优化步骤。

目标：

Vina
 ↓
protein side chain固定
 ↓
PyRosetta
 ↓
重新选择合理 rotamer

最基本：

task = pyrosetta.standard_packer_task(
    pose
)

task.restrict_to_repacking()

然后：

packer = (
    protocols.minimization_packing
    .PackRotamersMover(
        scorefxn,
        task
    )
)

packer.apply(
    pose
)

二十二、但是为什么不能全蛋白 repack？

因为：

Protein = 500 aa

如果：

500 aa全部repack

可能导致：

大量不必要的构象变化

所以科研中更合理：

Ligand
 ↓
6 Å
 ↓
允许 repack

而：

远离 ligand
 ↓
禁止 repack

二十三、Step 11：建立局部 TaskFactory

from pyrosetta.rosetta.core.pack.task import TaskFactory

from pyrosetta.rosetta.core.pack.task.operation import (
    InitializeFromCommandline,
    RestrictToRepacking,
    PreventRepackingRLT,
    OperateOnResidueSubset
)

建立：

tf = TaskFactory()

tf.push_back(
    InitializeFromCommandline()
)

tf.push_back(
    RestrictToRepacking()
)

然后：

tf.push_back(
    OperateOnResidueSubset(
        PreventRepackingRLT(),
        pocket_selector,
        flip_subset=True
    )
)

这里：

flip_subset=True

意味着：

pocket 之外禁止 repack。

二十四、创建 task

task = (
    tf.create_task_and_apply_taskoperations(
        pose
    )
)

然后：

packer = (
    protocols.minimization_packing
    .PackRotamersMover(
        scorefxn,
        task
    )
)

packer.apply(
    pose
)

这时候完成：

Side-chain optimization

二十五、Step 12：为什么还要 Minimize？

Repack 是：

离散 rotamer 搜索

而 Minimize 是：

连续坐标优化

也就是说：

Repack
 ↓
选择一个更好的 rotamer

Minimize
 ↓
在这个 rotamer 附近继续微调

所以：

Repack
 →
Minimize

是非常自然的组合。

二十六、Step 13：建立 MoveMap

movemap = MoveMap()

movemap.set_bb(False)

movemap.set_chi(False)

意思：

默认：

backbone = 不动
side chain = 不动

然后：

for resi in pocket_residues:

    movemap.set_chi(
        resi,
        True
    )

变成：

Protein backbone
      ↓
    固定

Pocket side chains
      ↓
    可以动

二十七、Step 14：Ligand 的自由度

对于 ligand：

movemap.set_chi(
    ligand_res,
    True
)

让 ligand 的可优化内部自由度参与优化。

但是这里要强调：

是否、以及如何让 ligand torsion 被 Rosetta 正确优化，取决于 .params 中 ligand 的 torsion 定义和具体 Rosetta/PyRosetta 版本。

因此不能简单认为：

set_chi(ligand_res, True)

就等于“完整 ligand docking”。

二十八、Step 15：MinMover

min_mover = (
    protocols.minimization_packing
    .MinMover()
)

min_mover.movemap(
    movemap
)

min_mover.score_function(
    scorefxn
)

min_mover.min_type(
    "lbfgs_armijo_nonmonotone"
)

min_mover.tolerance(
    0.01
)

然后：

min_mover.apply(
    pose
)

完成：

local energy minimization

二十九、这一阶段到底发生了什么？

假设：

Vina：

Phe
 \
  \
   LIG
  /
Asp

经过 Repack：

Phe ↻
 \
  \
   LIG
  /
Asp ↻

再经过 Minimize：

Phe
  \
   LIG
  /
Asp

把：

距离
角度
torsion
van der Waals
electrostatics

进一步调整到局部能量较低的状态。

三十、Step 16：FastRelax

接下来可以做局部 Relax。

但是：

FastRelax 是强力工具，不建议没有约束地对整个蛋白做。

如果只是想学习流程，可以先：

relax = protocols.relax.FastRelax()

relax.set_scorefxn(
    scorefxn
)

relax.apply(
    pose
)

但科研版本建议加入：

MoveMap
+
Coordinate constraints
+
Pocket restriction

避免蛋白整体漂移。

三十一、为什么需要 Coordinate Constraint？

例如 ligand 附近的蛋白：

Phe100
Tyr105
Asp110

允许它们：

side chain调整

但不希望：

Phe100整个跑掉

所以：

允许局部优化
+
限制整体偏移

这就是 coordinate constraint 的意义。

三十二、Step 17：再次 Minimize

Relax 后：

min_mover.apply(
    pose
)

再次进行局部最小化。

因此完整过程：

Vina
 ↓
Repack
 ↓
Minimize
 ↓
Relax
 ↓
Minimize

三十三：完整代码框架

import pyrosetta

from pyrosetta import *
from pyrosetta.rosetta import *

from pyrosetta.rosetta.core.select import residue_selector

from pyrosetta.rosetta.core.pack.task import (
    TaskFactory
)

from pyrosetta.rosetta.core.pack.task.operation import (
    InitializeFromCommandline,
    RestrictToRepacking,
    PreventRepackingRLT,
    OperateOnResidueSubset
)


# =====================================================
# 1. Init
# =====================================================

pyrosetta.init(
    "-extra_res_fa params/LIG.params "
    "-mute all"
)


# =====================================================
# 2. Load complex
# =====================================================

pose = pose_from_pdb(
    "complexes/complex_01.pdb"
)


# =====================================================
# 3. Find ligand
# =====================================================

ligand_res = None

for i in range(
    1,
    pose.total_residue() + 1
):

    if pose.residue(i).name3() == "LIG":

        ligand_res = i
        break


if ligand_res is None:

    raise RuntimeError(
        "LIG not found"
    )


print(
    "Ligand:",
    ligand_res
)


# =====================================================
# 4. Binding pocket
# =====================================================

lig_selector = (
    residue_selector.ResidueIndexSelector(
        str(ligand_res)
    )
)


pocket_selector = (
    residue_selector.NeighborhoodResidueSelector()
)

pocket_selector.set_focus(
    lig_selector
)

pocket_selector.set_distance(
    6.0
)

pocket_selector.set_include_focus_in_subset(
    True
)


subset = pocket_selector.apply(
    pose
)


pocket_residues = []

for i in range(
    1,
    pose.total_residue() + 1
):

    if subset[i]:

        pocket_residues.append(i)


print(
    "Pocket:",
    pocket_residues
)


# =====================================================
# 5. Score
# =====================================================

scorefxn = get_fa_scorefxn()

score0 = scorefxn(
    pose
)

print(
    "Initial:",
    score0
)


# =====================================================
# 6. Repack pocket
# =====================================================

tf = TaskFactory()

tf.push_back(
    InitializeFromCommandline()
)

tf.push_back(
    RestrictToRepacking()
)

tf.push_back(
    OperateOnResidueSubset(
        PreventRepackingRLT(),
        pocket_selector,
        flip_subset=True
    )
)


task = (
    tf.create_task_and_apply_taskoperations(
        pose
    )
)


packer = (
    protocols.minimization_packing
    .PackRotamersMover(
        scorefxn,
        task
    )
)

packer.apply(
    pose
)


score1 = scorefxn(
    pose
)

print(
    "After repack:",
    score1
)


# =====================================================
# 7. MoveMap
# =====================================================

movemap = MoveMap()

movemap.set_bb(
    False
)

movemap.set_chi(
    False
)


for resi in pocket_residues:

    movemap.set_chi(
        resi,
        True
    )


movemap.set_chi(
    ligand_res,
    True
)


# =====================================================
# 8. Minimization
# =====================================================

min_mover = (
    protocols.minimization_packing
    .MinMover()
)

min_mover.movemap(
    movemap
)

min_mover.score_function(
    scorefxn
)

min_mover.min_type(
    "lbfgs_armijo_nonmonotone"
)

min_mover.tolerance(
    0.01
)

min_mover.apply(
    pose
)


score2 = scorefxn(
    pose
)

print(
    "After minimize:",
    score2
)


# =====================================================
# 9. Relax
# =====================================================

relax = protocols.relax.FastRelax()

relax.set_scorefxn(
    scorefxn
)

relax.apply(
    pose
)


score3 = scorefxn(
    pose
)

print(
    "After relax:",
    score3
)


# =====================================================
# 10. Final minimize
# =====================================================

min_mover.apply(
    pose
)


score4 = scorefxn(
    pose
)

print(
    "Final:",
    score4
)


# =====================================================
# 11. Save
# =====================================================

pose.dump_pdb(
    "optimized/complex_01_refined.pdb"
)

三十四、但这还不是最终科研版

这份代码主要理解：

Pose
↓
ResidueSelector
↓
TaskFactory
↓
Repack
↓
MoveMap
↓
MinMover
↓
FastRelax

真正用于系统计算，建议再加入：

① 多个 Vina poses
Top 10
② 每个 pose 多次 refinement

例如：

10 Vina poses
×
10 Rosetta trajectories
=
100 structures
③ Pose clustering

例如：

100 structures
        ↓
RMSD clustering
        ↓
Cluster 1
Cluster 2
Cluster 3
④ Interface score

不要只看：

total_score

而要重点看：

interface score
interface ΔG

三十五、为什么需要多次 PyRosetta refinement？

因为 Rosetta optimization 也不是：

一次运行
↓
绝对全局最优

它仍然可能受到：

初始构象
rotamer
局部能量极小值

影响。

所以：

Vina Pose 1
    ↓
Rosetta 1
Rosetta 2
Rosetta 3
...

得到：

多个局部 minimum

如果最终很多结果都集中到：

相似 binding mode

那么这个 binding mode 更值得关注。

三十六、推荐的最终数据结构

比如：

Vina pose	Rosetta run	Vina score	Rosetta score	Interface ΔG	H-bond	Clash
P1	1	-10.2	-145	-25	5	0
P1	2	-10.2	-149	-27	6	0
P1	3	-10.2	-143	-24	4	0
P2	1	-9.9	-158	-31	7	0
P2	2	-9.9	-161	-33	8	0
P3	1	-9.6	-140	-20	3	1

会发现：

Vina最好
P1

不一定：

Rosetta refinement最好

可能：

P2

经过局部优化后 interface 更好。

三十七、最后一步：不要马上认为 Rosetta 最低分就是正确答案

正确的筛选应该是：

Vina score
       +
Rosetta interface energy
       +
H-bond
       +
vdW
       +
electrostatics
       +
steric clash
       +
pose clustering

最后再进入：

GROMACS MD

三十八、整个方法的科学逻辑

可以把这套方法理解成三个层次：

Level 1：Vina

解决：

Ligand 在哪里？

搜索 binding pose

↓

Level 2：PyRosetta

解决：

Protein-ligand interface 怎么调整更合理？

side chain
+
torsion
+
local geometry
+
interface energy

↓

Level 3：GROMACS

解决：

这个 complex 在动态环境中是否稳定？

protein dynamics
+
ligand dynamics
+
water
+
ions

这三个步骤实际上对应三个不同的问题：

Vina：能不能放进去？

PyRosetta：放进去以后，界面能不能调整得更合理？

GROMACS：调整后的复合物在动态环境中能不能稳定存在？

最后提醒一点：PyRosetta 局部优化并不能保证把错误的 Vina pose “纠正成正确 pose”。如果 Vina 一开始把 ligand 放到了完全错误的 binding mode，local refinement 通常只能在这个局部构象附近优化。因此科研上更推荐 Top 5–10 Vina poses → 每个 pose 多次 PyRosetta refinement → clustering + interface 分析 → 再选 Top 2–3 进入 MD，而不是只优化 Vina 的第一名。

如果接下来要实际运行，建议下一步直接解决最关键的技术环节：写一套完整的 Vina PDBQT → 自动拆分 Top 10 pose → 保留 Vina 坐标 → 自动生成/匹配 LIG.params → 自动构建 protein-ligand complex → PyRosetta 局部 Repack/Minimize/FastRelax → 输出 CSV 排名的可运行 pipeline。这样才是真正从 Vina 结果一路跑到 PyRosetta，而不是停留在示例代码层面。

生活很好，有你更好

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

数据分析

#分子对接

#分子动力学

#数据分析

作者，Evil Genius

接下来要为我们的生化课程，分子对接/分子动力学做好准备，大家学习分子对接/分子动力学，autodock vina、pyrosetta、Gromacs是必要要掌握的。

这一篇我们来采用pyrosetta对autodock vina的分子对接结果进行优化/

首先了解一下背景
- Vina 不是简单的“刚性对接”。
  - 更准确地说，经典 AutoDock Vina 是 rigid receptor + flexible ligand：受体通常保持固定，而 ligand 的平移、旋转和可旋转键构象可以搜索。因此，Vina 最大的局限之一不是“配体不能动”，而是蛋白结合口袋的柔性没有被充分考虑。Vina 官方文档也将 receptor rigidity 与 ligand flexibility 分开处理。
- 一、首先理解：Vina 到底在做什么？
  - 假设蛋白：
  - ligand：
- Vina 要解决的问题是：
  - 在指定 binding box 中，寻找 ligand 比较合理的空间位置和构象。
  - 可以理解为：
  - 最后得到：

二、Vina 的第一个核心局限：受体通常是刚性的
- 假设真实结合过程：
  - 未结合：
  - Ligand进入：
  - 真实情况下：
  - ligand 进入 binding pocket 后，蛋白侧链可能发生 rotamer rearrangement。
- 但是经典 Vina：
  - 也就是说：
  - 所以它更准确的定义是：
  - rigid-receptor / flexible-ligand docking
- 而不是：
  - rigid ligand docking。

三、第二个局限：Vina 的 binding pocket 不会充分“诱导适配”
- 蛋白-配体结合经常存在：
  - 例如：
  - 也就是：
- Vina 的普通 rigid receptor 模式不能充分描述这一过程。
- 因此：
  - Vina 得到的是一个 docking pose，而不一定是最终能量最低的 protein-ligand complex。

四、第三个局限：Vina 的最佳 pose 不一定是真实 pose
- 例如：
- 不能直接认为：
  - 因为 docking score 是一个近似评分函数。
  - 它并不是严格的：
  - 实验结合自由能
- 而是：
  - 近似 scoring function
  - 所以可能出现：
- 而：
- 因此：
  - 不能只看 Vina score。

五、第四个局限：蛋白侧链构象可能来自 apo 状态
- 比如蛋白 PDB 是：
- 那么：
  - 的构象可能是：
- 但是 ligand 结合后：
  - 可能不同。
  - 于是：
- 得到的结果可能受到：
  - 初始蛋白构象的限制。

六、第五个局限：Vina 不是 MD
- Vina 做的是：
- 而不是：
- 所以它不会告诉：
- 这些问题应该由：
  - 进一步解决。

七、那么为什么使用 PyRosetta？
- 这里就出现了一个非常自然的逻辑：
- 所以：
  - Vina 负责“找到一个合理的结合姿势”，PyRosetta 负责“在这个姿势附近把蛋白-配体界面调整得更合理”。

八、为什么一定强调“局部优化”？
- 因为并不是想重新折叠整个蛋白。
- 假设：
- ligand 只和：
  - 发生直接相互作用。
- 那么真正关心的是：
- 而不是：
- 因此：
- 这就是：

九、完整工作流程
- 最终推荐：

十、Step 1：准备文件
- 建议目录：

十一、Step 2：Vina docking
- 例如：
- 这里：
- 让 Vina 输出多个 pose。
  - 建议：
  - 作为最低限度。
- 更稳妥：

十二、Step 3：为什么保留 Top 10？
- 因为：
- 不是绝对真值。
  - 例如：
- 这些 pose 可能代表：
  - 同一个 binding mode 的不同构象
- 也可能代表：
  - 不同 binding modes
- 所以后续 PyRosetta 应该：
- 而不是：
  - 只 refine Pose 1

十三、Step 4：准备 ligand 参数
- 这是 PyRosetta 最重要的准备工作之一。
- Vina 使用：
- PyRosetta 使用：
- 所以：
- Rosetta 官方的 ligand preparation 流程就是通过 molfile_to_params.py 从分子文件生成 ligand params。
  - 例如：
  - 得到：

十四、.params 到底是什么？
- 可以把它理解成：
  - 告诉 Rosetta：“LIG 是什么东西”。
- 包括：
  - 因此：
- 只有坐标是不够的。
- Rosetta 还需要：
  - 来理解 ligand 的化学结构。

十五、Step 5：Vina pose 和 Rosetta ligand 坐标怎么结合？
- 这里是整个 pipeline 中最容易出错的部分。
- 实际上有两个信息来源：
- 化学拓扑
  - 来自：
- ↓
- 空间坐标
  - 来自：
- 所以：
- 最终：
- 组成：

不要简单地把 PDBQT 改成 PDB 后就结束。

十六、Step 6：PyRosetta 初始化
- 然后：
- 加载：
- 检查：
- 应该看到：

十七、Step 7：找到 ligand

十八、Step 8：找 binding pocket
- 假设使用：
  - 作为局部优化范围。
- 得到 pocket：

十九、为什么是 6 Å？
- 不是一个绝对规则。
  - 可以尝试：
  - 通常：
  - 更严格。
  - 是一个比较实用的起点。
  - 允许更大范围的局部适配。
- 对于第一次建立 pipeline：
  - 推荐先使用 6 Å。

二十、Step 9：计算初始 Rosetta score
- 注意：
  - 这个 score 和 Vina score 不能直接比较。
  - 例如：
- 不能说：
- -135 比 -10.2 好
  - 它们是完全不同的 scoring function。

二十一、Step 10：Side-chain Repacking
- 这是第一个真正的优化步骤。
  - 目标：
  - 最基本：
  - 然后：

二十二、但是为什么不能全蛋白 repack？
- 因为：
  - 如果：
  - 500 aa全部repack
- 可能导致：
  - 大量不必要的构象变化
- 所以科研中更合理：
- 而：

二十三、Step 11：建立局部 TaskFactory
- 建立：
- 然后：
- 这里：
- 意味着：
  - pocket 之外禁止 repack。

二十四、创建 task
- 然后：
  - 这时候完成：

二十五、Step 12：为什么还要 Minimize？
- Repack 是：
  - 离散 rotamer 搜索
- 而 Minimize 是：
  - 连续坐标优化
- 也就是说：
- 所以：
  - 是非常自然的组合。

二十六、Step 13：建立 MoveMap
- 意思：
- 默认：
- 然后：
- 变成：

二十七、Step 14：Ligand 的自由度
- 对于 ligand：
  - 让 ligand 的可优化内部自由度参与优化。
- 但是这里要强调：
  - 是否、以及如何让 ligand torsion 被 Rosetta 正确优化，取决于 .params 中 ligand 的 torsion 定义和具体 Rosetta/PyRosetta 版本。
- 因此不能简单认为：
  - 就等于“完整 ligand docking”。

二十八、Step 15：MinMover
- 然后：
- 完成：

二十九、这一阶段到底发生了什么？
- 假设：
  - 经过 Repack：
  - 再经过 Minimize：
- 把：
  - 进一步调整到局部能量较低的状态。

三十、Step 16：FastRelax
- 接下来可以做局部 Relax。
- 但是：
  - FastRelax 是强力工具，不建议没有约束地对整个蛋白做。
  - 如果只是想学习流程，可以先：
- 但科研版本建议加入：
  - 避免蛋白整体漂移。

三十一、为什么需要 Coordinate Constraint？
- 例如 ligand 附近的蛋白：
  - 允许它们：
- 但不希望：
  - 所以：
  - 这就是 coordinate constraint 的意义。

三十二、Step 17：再次 Minimize
- Relax 后：
- 再次进行局部最小化。
- 因此完整过程：

三十三：完整代码框架

三十四、但这还不是最终科研版
- 这份代码主要理解：
- 真正用于系统计算，建议再加入：
- 例如：
- 例如：
- 不要只看：
- 而要重点看：

三十五、为什么需要多次 PyRosetta refinement？
- 因为 Rosetta optimization 也不是：
- 它仍然可能受到：
- 影响。
- 所以：
- 得到：
- 如果最终很多结果都集中到：
- 那么这个 binding mode 更值得关注。

三十六、推荐的最终数据结构
- 比如：

会发现：
- 不一定：
- 可能：
  - 经过局部优化后 interface 更好。

三十七、最后一步：不要马上认为 Rosetta 最低分就是正确答案
- 正确的筛选应该是：
  - 最后再进入：

三十八、整个方法的科学逻辑
- 可以把这套方法理解成三个层次：
- Level 1：Vina
  - 解决：
  - Ligand 在哪里？
- ↓
- Level 2：PyRosetta
  - 解决：
  - Protein-ligand interface 怎么调整更合理？
- ↓
- Level 3：GROMACS
  - 解决：
  - 这个 complex 在动态环境中是否稳定？

这三个步骤实际上对应三个不同的问题：
- Vina：能不能放进去？
- PyRosetta：放进去以后，界面能不能调整得更合理？
- GROMACS：调整后的复合物在动态环境中能不能稳定存在？
- 最后提醒一点：PyRosetta 局部优化并不能保证把错误的 Vina pose “纠正成正确 pose”。如果 Vina 一开始把 ligand 放到了完全错误的 binding mode，local refinement 通常只能在这个局部构象附近优化。因此科研上更推荐 Top 5–10 Vina poses → 每个 pose 多次 PyRosetta refinement → clustering + interface 分析 → 再选 Top 2–3 进入 MD，而不是只优化 Vina 的第一名。
- 如果接下来要实际运行，建议下一步直接解决最关键的技术环节：写一套完整的 Vina PDBQT → 自动拆分 Top 10 pose → 保留 Vina 坐标 → 自动生成/匹配 LIG.params → 自动构建 protein-ligand complex → PyRosetta 局部 Repack/Minimize/FastRelax → 输出 CSV 排名的可运行 pipeline。这样才是真正从 Vina 结果一路跑到 PyRosetta，而不是停留在示例代码层面。

生活很好，有你更好