第3章 基础输入文件：PDBQT格式详解、受体和配体的预处理
做分子对接，输入文件就是你的“食材”。食材不对，厨艺再好也白搭。
AutoDock Vina用的不是普通的PDB文件，而是PDBQT格式。这个格式多了一个Q和T。Q代表电荷（partial charge），T代表可旋转键类型（torsion type）。说白了，它就是给普通PDB文件加上了分子力学必需的电性信息和柔性信息。
3.1 PDBQT格式的四个核心模块
打开一个PDBQT文件，你会发现它分成了四个部分。我刚开始学的时候，觉得这格式挺啰嗦的。后来才明白，每一行都有它存在的意义。
| 模块 | 作用 | 典型行数 |
|---|---|---|
| REMARK | 文件元信息，生成工具、版本号、来源PDB ID | 5-20行 |
| ROOT | 刚性部分的原子坐标和参数，不能旋转 | 按原子数定 |
| BRANCH | 可旋转分支上的原子，每个分支单独标记 | 按分支数定 |
| TORSDOF | 末尾一行，记录总的可旋转键自由度 | 1行 |
看看下面这个配体PDBQT的片段，就更容易理解了：
REMARK Name = ligand
REMARK Active torsions: 5
ROOT
ATOM 1 C1 LIG 1 1.234 2.345 3.456 0.150 0.000
ATOM 2 C2 LIG 1 1.567 2.890 3.123 0.100 0.000
...
ENDROOT
BRANCH 3 4
ATOM 4 C4 LIG 1 1.901 3.234 2.789 0.050 0.000
...
ENDBRANCH
TORSDOF 5每一行原子的格式是：ATOM序号 原子名 残基名 链ID 残基序号 X Y Z 电荷 类型
最后一个字段是AutoDock专用的原子类型，不是元素符号。比如C是芳香碳，A是极性碳，N是氮，OA是氧等等。
原子类型速查表（我常用的几个）：
- C – 非极性碳（烷基、芳香碳骨架）
- A – 极性碳（连有杂原子的碳）
- N – 氮
- OA – 氧（羰基、羟基中的氧）
- HD – 可形成氢键的氢
- P – 磷
- S – 硫
- F – 氟、Cl – 氯、Br – 溴、I – 碘
3.2 受体预处理：刚体到底准备什么
受体在对接时默认是刚性的，也就是说所有原子固定不动。但这不代表你随便拿个PDB就能用。
我在项目中遇到过好几次，从PDB库里下载的受体结构，缺失了好几个关键残基的侧链。直接去对接，结果可想而知，活性位点都不完整。
受体预处理的步骤，我来列个清单：
- 去除水分子：除非你知道某个水分子参与了关键的氢键网络，否则全部删掉。
- 补齐缺失残基/原子：用Modeller、PDBFixer等工具补全。这一步不能省。
- 加氢：调整pH 7.4下的质子化状态。His残基的侧链常见3种状态，选错了会直接影响对接结果。
- 计算电荷：分配Gasteiger电荷，非极性氢合并。
- 输出PDBQT：用MGLTools或ADFR的脚本生成。
来看看我推荐用的命令：
# 使用 prepare_receptor4.py (MGLTools)
/opt/mgltools/bin/pythonsh /opt/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py \
-r receptor.pdb \
-o receptor.pdbqt \
-A hydrogens \
-U nphs_lps_waters
# -A hydrogens：加氢
# -U nphs_lps_waters：去除非极性氢、孤对电子、水有个小技巧：先打开PDB文件看一下，有些PDB里已经自带氢原子了，这时候再用-A hydrogens可能会造成多余加氢。我习惯先检查REMARK里有没有“HETATM”水分子以外的氢原子记录。有的话就用 -A None。
3.3 配体预处理：电荷与可旋转键的坑
配体比受体麻烦得多。因为你得决定哪些键可以转，哪些不能转。这个决策直接影响对接的效率和质量。
配体处理的标准流程：
- 从ChemDraw、PubChem等来源得到SDF或MOL2格式
- 用Open Babel转换为PDB（或直接MOL2）
- 用prepare_ligand4.py分配原子类型、计算电荷、检测可旋转键
- 输出PDBQT
# 使用 prepare_ligand4.py
/opt/mgltools/bin/pythonsh /opt/mgltools/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_ligand4.py \
-l ligand.pdb \
-o ligand.pdbqt \
-F \
-v
# -F：允许检测单键作为可旋转键
# -v：verbose模式，看脚本在做什么电荷设置，这里水很深。prepare_ligand4.py默认分配Gasteiger电荷。这个电荷模型速度快，对中性小分子问题不大。但碰到带电配体，比如季铵盐、磷酸基团，Gasteiger电荷可能会失真。
我曾经就吃过大亏。一个带两个磷酸基团的配体，用Gasteiger电荷跑完后，对接能量怎么都不对。后来换成基于DFT的RESP电荷，结果就正常了。
# 手动指定配体电荷（假设你已经有了mol2格式带RESP电荷的文件）
# 可以直接用如下命令生成PDBQT，跳过电荷重新计算
obabel ligand_charged.mol2 -O ligand.pdb
# 然后用prepare_ligand4.py处理，-U charge参数可以跳过电荷计算
prepare_ligand4.py -l ligand.pdb -o ligand.pdbqt -U charge避坑指南：
- 对于含金属离子的配体，AutoDock的原子类型可能不识别，需要手动查阅ADT原子类型列表，选择最近的类型。
- 环内单键尽量不要设为可旋转键。程序有时会把环上的单键也识别为可旋转，你要手动改“TORSDOF”行或调整-t参数。
- 肽类配体（多肽）的处理很特殊。建议直接用MGLTools的GUI看哪个键该设成可旋转，不要全自动。
3.4 可旋转键的检测背后的逻辑
说白了，程序就是在数配体里的单键数量。但不是所有单键都能转。双键不能转，三键不能转，环里的单键不能转，末端甲基的C-H键也不用转。
你想想看，prepare_ligand4.py怎么判断一个键该不该保留为刚性？它用的是一套“原子杂化类型+连接拓扑”的规则集。我在做批量对接时，经常发现同一个系列的化合物，脚本判断出的可旋转键数却不一样。原因往往是某个取代基上的sp3碳连到了sp2碳上，导致键的类型判断漂移。
下面是可旋转键判断的简化逻辑，用表格展示：
| 键类型 | 是否可旋转 | 说明 |
|---|---|---|
| C-C单键（两端都是sp3） | 是 | 典型的可旋转键 |
| C-C单键（一端sp3，一端sp2） | 是 | 但旋转会改变共轭，慎用 |
| C=C双键 | 否 | 刚性，不可旋转 |
| C≡C三键 | 否 | 刚性 |
| 环内单键 | 否 | 刚性 |
| C-N单键（酰胺） | 否 | 酰胺有部分双键特征，刚性 |
| O-H、N-H | 否 | 末端键，不管 |
3.5 知识结构流程图
我把PDBQT生成的核心流程画成了下面这张SVG图，你看一遍就清楚了：
3.6 批量处理时参数标准化
很多人用Vina做批量对接时，预处理这一步做得很随意。有的受体加了氢，有的没加；有的配体电荷计算错了；有的可旋转键数目忽高忽低。这样的结果放在一起，能比吗？
批量预处理参数标准化的核心原则：
- 统一加氢工具和版本（比如全部用MGLTools 1.5.7）
- 统一电荷模型（Gasteiger还是RESP，统一选一个）
- 统一去除溶剂策略（结合水保留还是全删？）
- 统一可旋转键数目的上限（比如最多不超过10个）
- 统一输出命名规则（比如 受体_chain.pdbqt, 配体_ID.pdbqt）
我曾经在做一个包含2000个化合物的虚拟筛选项目时，一开始没注意配体的可旋转键数目标准化。跑完之后发现，有些结果太好，我都不敢信。后来一查，原来是那些配体被错误地设成了全刚性，根本不科学。排查了两天才找到原因。
从那以后，我都在预处理流程里加一行检查命令：
# 批量检查配体PDBQT中的可旋转键数
grep "TORSDOF" *.pdbqt | awk '{print $1, $2}'跑完一看，所有配体都能看到它们的自由度。超过10的我会单独拿出来，看看是不是有不该转的键被误判了。这一个简单的grep命令，帮我省了太多事后排查的麻烦。
标准化检查清单：
- 受体PDBQT是否包含完整活性位点？（缺失残基？）
- 配体PDBQT的原子类型是否正确？（尤其注意金属、卤素）
- 可旋转键数是否在合理范围内？（小分子一般3-8个）
- 电荷总和是否为整数？（配体电荷应为0或±1/±2等整数）
- 文件命名是否统一？（推荐：受体名称_配体名称.pdbqt）
嗯，预处理做到位了，你的对接结果才靠谱。这就是所谓的“垃圾进，垃圾出”。Garbage In, Garbage Out. 没商量。