SPONGE常用教程：蛋白+配体模拟1

软件支持
SPONGE（Simulation Package tOward Next GEneration molecular modelling）是由北京大学高毅勤课题组开发的分子动力学模拟程序。安装教程
XPONGE 使用python编写的分子动力学模拟前后处理工具。简易安装：pip install xponge
DSDP 使用GPU开发的蛋白-分子对接工具。

目录

1. 应用场景简述；
2. DSDP:蛋白-配体对接；
3. XPONGE：蛋白-配体建模，加溶剂；
4. SPONGE：能量极小化-NVT-NPT-正式模拟；
5. XPONGE：数据简单后处理。

1. 应用场景简述
研究蛋白与药物分子结合的体系中，模拟盒子包含以下几个部分：A.蛋白 + B.有机配体分子 + C. 溶剂水分子。这一过程包括以下几个步骤：
a. DSDP:蛋白-配体对接(可选)如果蛋白-配体之间不在合适的位置，可以使用分子对接算法将两者坐标调整到正确的位置上。输入：蛋白pdb+配体mol2；输出：蛋白pdb+配体mol2（对接坐标）
b. XPONGE:蛋白-配体建模，加溶剂（必选）将蛋白、配体和溶剂分别加载对应力场，加入到体系中，生成SPONGE模拟引擎可读文件。输入：蛋白pdb+配体mol2（对接坐标）；输出：一系列文本文件+体系pdb
c. SPONGE:分子动力学模拟（必选）体系经过能量极小化优化氢原子位置，升温到反应温度，加压到反应压强，然后运行成品模拟。输入：一系列文本文件；输出：轨迹（二进制文件）
d. XPONGE:数据简单后处理（可选）统计轨迹中的体系稳定性信息（RMSD、Rg、RMSF）及氢键数等基础信息。输入：轨迹+体系pdb；输出：统计文件txt

数据预清洗(1AE5-HUP作用对,点击下载原文件)：
PDB数据库中下载的PDB文件（1ae5_pdbraw.pdb）除了蛋白的PDB结构之外，往往还会带有额外稳定结构的配体或者结晶水等我们不需要的信息，需要进行预清洗。
这里我们使用Xponge.pdb_filter函数来获得简化的文件，其中前两个参数分别是输入和输出的文件名，heads指定需要的表头部分，hetero_residues指定需要的非蛋白残基（这个例子没有）。

import Xponge
Xponge.pdb_filter("1ae5_pdbraw.pdb", "1ae5_protein.pdb", heads=["ATOM", "SEQRES", "TER"], hetero_residues=[])

小分子配体mol2文件(HUP.mol2)有可能会存在氢原子缺失的情况，若是，可使用openbabel对mol2进行加氢处理，设置pH为7。

obabel 1ae5_HUP_ligand.mol2 -O HUP_ligand.mol2 -p 7

确认蛋白-配体相对位置：
可通过分子可视化软件（例如pymol,vmd等）确认蛋白与配体的相对位置，坐标是否匹配，分以下两种情况：

如图左，绿色为蛋白6OOY，紫色为配体A7M(C17H18N2O)的正确结合位置，可以跳过第2步，直接进行第3步； 图右，蛋白1AE5与配体HUP(紫色，Huperzine A，C15H18N2O)有相互作用，但是未找到对接的pdb结构，需要执行DSDP分子对接。

2. DSDP:蛋白-配体对接

下面以1AE5-HUP作用对为例，演示对接操作：
安装DSDP
**【安装openbabel格式转换工具】（pip install openbabel，第三方工具，导出格式不能保证正确可读，建议用可视化软件打开进行double-check）
**【安装MGLtool格式转换工具】(参考https://ccsb.scripps.edu/mgltools/downloads/,第三方工具，用于生成pdbqt文件，常出现键型信息错误，同样需要校验)

conda activate DSDP
#DSDP成功安装后，进入DSDP_python环境
cd DSDP
#进入DSDP文件夹，复制MGLtools可执行文件pythonsh、转换脚本prepare_receptor4.py和prepare_ligand4.py到DSDP所在的文件夹中进行操作：

mkdir run_pdb
cd run_pdb
mkdir 1ae5
#根据蛋白名新建文件夹
cd 1ae5
cp ../../../1ae5_protein.pdb .
cp ../../../HUP_ligand.mol2 .
#将蛋白文件（1ae5_protein.pdb）和配体文件（HUP_ligand.mol2）复制到当前文件夹
../../pythonsh ../../prepare_receptor4.py -r 1ae5_protein.pdb -A hydrogens -o 1ae5_protein.pdbqt
../../pythonsh ../../prepare_ligand4.py -l HUP_ligand.mol2 -A hydrogens -o 1ae5_ligand.pdbqt
#此时文件夹中会多出1ae5_protein.pdbqt和1ae5_ligand.pdbqt两个文件
cd ../../

MGLtools为第三方工具，安装过程中可能出现报错，或者找不到正确文件等问题。这里提供生成的pdbqt文件，复制到1ae5文件夹中，直接进入下一步。

#位于DSDP文件夹
python DSDP_blind_docking.py \
--dataset_path ./run_pdb/ \
--dataset_name run_pdb \
--site_path ./results/DSDP_dataset/site_output/ \
--exhaustiveness 384 \
--search_depth 40 \
--top_n 1 \
--out ./results/DSDP_dataset/docking_results/ \
--log ./results/DSDP_dataset/docking_results/
#定位到生成的文件夹中
cd results/DSDP_dataset/docking_results/
#此处需要通过其他软件转换pdbqt文件回到mol2格式，教程采用openbabel进行格式转换
obabel -ipdbqt 1ae5_out.pdbqt -opdb -O 1ae5_out.pdb
#查看配体pdb结构，发现-NH2基团键级错误（来自autodock-MGLtool的prepare_ligand4.py脚本在格式转换过程中误识别成-NH3+），补救措施为手动删除1ae5_out.pdb的id=22的氢原子(删除第51&29行，再删去第47行的“22  ”字符，且第40行CONNECT应该多加一个“13”，第41行CONNECT应该多加一个“12”，若有更优的解决方法欢迎讨论)，再转mol2格式
obabel -ipdb 1ae5_out.pdb -omol2 -O 1ae5_HUP_ligand.mol2 -p 10
#通常加氢应该是加到ph=7，此处配体分子式（C15H18N2O）对应-NH2，而非-NH3+，故采用ph=10保存配体分子
cp 1ae5_HUP_ligand.mol2 ../../../../
#将pdbqt文件转回mol2文件，复制到工作文件夹

再次查看配体位置，发现已经对接到蛋白口袋。

3.XPONGE：蛋白-配体建模，加溶剂
将蛋白文件（1ae5_protein.pdb）和配体文件（1ae5_HUP_ligand.mol2）复制到工作文件夹。【文件直接下载】

mkdir 1ae5
#创建空文件夹用于保存SPONGE输入文件（因为子文件数量较多（10+），放在一个文件夹目录下方便一起拷贝）

import Xponge
import Xponge.forcefield.amber.gaff as gaff
assign1 = Xponge.Get_Assignment_From_Mol2("1ae5_HUP_ligand.mol2")
assign1.Determine_Atom_Type("gaff")
assign1.Calculate_Charge("resp")
HUP = assign1.toResidueType("HUP")
gaff.parmchk2_gaff(HUP, "HUP.frcmod")
Xponge.save_mol2(HUP, "HUP.mol2")
#保存HUP分子力场用于后续调用

import Xponge.forcefield.amber.ff14sb
C = Xponge.load_pdb("1ae5_protein.pdb", ignore_hydrogen=False, ignore_unknown_name=True, ignore_seqres=True)
C.add_missing_atoms()
HUP = Xponge.load_mol2("HUP.mol2")
gaff.parmchk2_gaff(HUP, "HUP.frcmod")
a = C | HUP

#如果是真空体系的话，可以采用下面命令行保存，同样需要在外部创建空文件夹1ae5V
#Xponge.save_pdb(a, "1ae5_pro_lig.pdb")
#Xponge.save_sponge_input(a, "1ae5V/1ae5")

#正常模拟生理状态加水的话，用下面的命令行
import Xponge.forcefield.amber.tip3p
Xponge.addSolventBox(a, WAT, 5)
#此处加溶剂水为tip3p水，参数5代表蛋白距离盒子边界5A，通常模拟会选更大的距离，以避免边界效应
print(C.charge+int(round(HUP.charge)))
#此处显示体系带3个单位负电荷，需要加入离子补到电中性#思考题，若要达到生理盐水的离子浓度，此处电荷数量是怎样算出来的？
Xponge.Solvent_Replace(a, WAT, {CL:8})
Xponge.save_pdb(a, "1ae5_pro_lig_sol.pdb")
Xponge.save_sponge_input(a, "1ae5/1ae5")

	AtomId	ResId	ResCount
Protein	0-3351	0-222	223
Ligand	3352-3387	223	1
Iron	3388-3395	224-231	8
Water	3396-12650	232-3316	3085

[To be continued...]