Autodock vina の pdbqt准备

这篇blog中假设你毫无autodock经验但是有最最基础linux使用经验，本blog着重描述vina的一些文件和操作流程，以避免以后自己操作的时候犯错。

我们会先从配体准备的软件选择和安装开始，然后会分别讲解ligand和receptor的准备过程。

pdbqt生成-软件的选择和安装

这里的软件指的是从sdf或者pdb文件到pdbqt文件的过程。本质上来说他们都是笛卡尔坐标系xyz的信息储存文件，都可以用文本文档直接打开。

head your_pdbqt.pdbqt

如果你打开看了，就会发现其实pdbqt比pdb多出了两列：也就是q和t。Partial Charge (Q)和Atom Type (T)。Q用于计算静电相互作用，T用于计算van derr waal力和氢键，他们正常pdb是不需要的，但是预测docking的时候非常有用。

一直以来Autodock tools或者MGL tools是比较被推荐的，也是官方标准工具，也是图形化工具。另外也有PyRx等别的第三方pdbqt生成器的库，但是兼容性相对差一点点（生成QT就说白了是很简单的活儿，但是既然有官方的为啥不用呢）。但是对于跑在服务器上面，我还是更倾向于完全使用命令行完成，这样不管是文件的管理还是流程把控都比较容易，写sh脚本可以省去很多人工体力劳动还有无聊的scp和文件整理环节。这里强烈推荐meeko库（https://github.com/forlilab/Meeko）。这是autodock官方主力推动的，新的处理工具，对于vina兼容性是最好的。

利用pip install或conda install安装即可，你会获得mk_prepare_ligand.py 和mk_prepare_receptor.py 两个可以直接执行的脚本文件，安装完成以后可能还要手动安装一些第三方库如gemmi。整个过程比较简单可以和AI老师交流一下，由于大家的环境各不相同，我就不做详细说明了。

Ligands 准备

这里我们用类固醇Kaempferol来打个比方。一般来说都是chemdraw直接画好结构导出成sdf，然后进行处理。可以使用下面链接下载ligand的sdf

1 平衡pH (加H）

第一步是平衡到实验pH，第一步是平衡到实验pH，第一步是平衡到实验pH。曾经多次实验失败全都是因为没有注意pH。学化学的知道，H这玩意真是难管又重要。一定要记住这一步！！！所以第一步我们平衡pH

obabel Kaempferol.sdf -O Kaempferol_ph7.0.sdf -p 7.0

1.5 energy minimization

obabel Luteolin.sdf -O Luteolin_minimized.sdf --gen3d --minimize --ff MMFF94 --steps 500

有时候能量最小化一下，要不然可能不是真实情况。autodock会根据和输入的pdbqt的区别来扣分(如cis-trans转换等等)。我们需要模拟真实情况下的输入才能保证输出的准确性。

2 用meeko官方工具prepare ligand

下面直接用这个ph7.0的去prepare ligand即可，无需其他步骤。meeko会进行一些简单的基于经验的计算，正儿八经的算的话不会这么快的。

mk_prepare_ligand.py -i Kaempferol_ph7.0.sdf -o Kaempferol.pdbqt

这样我们就得到了Kaempferol.pdbqt!

receptor

autodock对于receptor的要求可谓是非常的严格，严格在于H，没错，又是H。为了降低工作量，vina只保留极性H，非极性H应该完全去除。另外，pH也是我们需要注意的，receptor使用3NPO做个例子。下载链接：https://www.rcsb.org/structure/3NPO

1 去水

这有很多种方法可以选，比如说用pymol打开手动选择所有0去掉，不过这里推荐直接一行命令搞定，不过还是要提前看一下水是不是HOH，如果是不讲科学不明来源的pdb，那这个方法可能不行。

grep -v HOH 3npo.pdb > 3npo_nowater.pdb

2 平衡pH，加H

途径1：reduce2

这里最为推荐的是autodock官方自家推出的reduce2.py，直接执行命令即可：

reduce2.py approach=add add_flip_movers=True 3npo_nowater.pdb > 3npo_charged.pdb

但是！！！！但是！！！他不支持滴定加H！！也就是说你不知道它加出来到底pH是多少！虽然来说一般都是ph 7~7.4，当然你也可以加完以后用propka进行预测，根据电荷判断，但是仍然无法保证蛋白就是你想要的质子化状态。

途径2：pdb2pqr

这里介绍一下另一个命令行工具pdb2pqr，这也是地位极高、神一样存在的脚本。他会把pdb文件变成pqr。pqr其实是pdbqr，也是在pdb的基础上加东西，这次加的两列是charge (Q)和radius (R)。

pdb2pqr30 --ff=AMBER --with-ph=7.0 3npo_nowater.pdb 3npo_ph7.0.pqr

这个命令行工具非常厉害，采用了propka算法，滴定计算pH环境，完全不是盲目用静态规则猜的。随后还会快速的进行几何优化，保证蛋白质H-bond网络总能量最低。但是对于autodock来说有点大材小用，如果你不需要确保ph或者reduce本身ph没问题的话。

如果你像我一样使用的是pdb2pqr，那么非常可惜非常可惜，pdb2pqr提供了q的信息并且还有算出来的r的信息，但是meeko不支持它作为输入，我们只能先丢掉这些信息 – 但是此时H已经被按照pH加上并且结构进行了一定程度的优化，所以也没白做。我们要把它转化成pdb再进行下一步

obabel 3npo_ph7.0.pqr -O 3npo_charged.pdb

3 meeko prepare receptor

到这一步，不管你用的什么途径，都会得到3npo_charged.pdb. 下面直接使用mk prepare即可

mk_prepare_receptor.py --read_pdb 3npo_charged.pdb -o receptor -p

有时候情况比较特殊我们需要注意一下二硫键, 整体的sh如下：

# 去掉水分子！
prody fetch 3npo
prody select "chain A and not water and not hetero" 3npo.pdb -o 3NPO_A.pdb
# ================= 如果有构象 ====================
# 1. 剔除所有带有 'B' 构象前缀的残基行 (如 BMET, BLEU, BGLU, BILE, BSER, BLYS)
grep -v -E "BMET|BLEU|BGLU|BILE|BSER|BLYS" 3NPO_A.pdb > 3NPO_cleaned.pdb
# 2. 将所有 'A' 构象前缀的残基名字改回标准的三字母缩写 (如 AMET ->  MET)
# 注意：前缀带有空格，用来保持 PDB 的严格列宽对齐
sed -i -E 's/AMET/ MET/g; s/ALEU/ LEU/g; s/AGLU/ GLU/g; s/AILE/ ILE/g; s/ASER/ SER/g; s/ALYS/ LYS/g' 3NPO_cleaned.pdb
# ========================================================
# mk prepare啦！确保这里double sulfur bond的template！
mk_prepare_receptor.py -i 3NPO_cleaned.pdb -o 3NPO -p -j -g --box_enveloping 3NPO_cleaned.pdb --padding 1 --set_template A:66,106,119,160=CYX
# 如果加柔性，加上下面的
# --flexres "A_ASN109,A_ASN90,A_MET107,A_LEU39,A_SER116"

这样我们就得到了receptor.pdbqt用于后续分析。

下一步？

你都拿到两个pdbqt啦！！直接写config file用vina dock就行！