1.5. AMBER案例#
Amber入门案例——蛋白质的折叠
教程来自 Amber 官网 basic tutorial3: Case Study - Folding TRP Cage (Advanced analysis and clustering), 展示了论文 "稳定蛋白折叠色氨酸肽笼的全原子结构预测和折叠模拟"中的工作.
这里以含 20 个氨基酸的肽 “trpcage” 为例, 模拟肽链的折叠过程.
1.5.1. 环境准备#
ambertools 工具的安装参见 安装GROMACS 节
pymol 软件用来查看肽链的结构
本节教程使用的所有文件打包在此: 百度网盘
1.5.2. 生成肽链#
创建一个输入文件 trpcage.in，用于生成初始的肽链, 包含以下内容：
# 加载力场和溶剂参数
source leaprc.protein.ff14SBonlysc # 加载ff14SBonlysc力场
set default PBradii mbondi3 # 设置GB-Neck2(隐式溶剂模型)所需的原子半径
# 构建肽链序列
trpcage = sequence { NASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO CSER }
# 检查结构完整性
Check trpcage
# 保存拓扑和坐标文件
savepdb trpcage trpcage_initial.pdb
saveamberparm trpcage trpcage.parm7 trpcage.crd
quit
然后使用tleap程序运行：tleap -f chignolin.in ,应该得到 tleap 消息:
...
Checking 'trp'....
/home/aiwt/miniconda3/bin/teLeap: Warning!
The unperturbed charge of the unit (1.000000) is not zero.
/home/aiwt/miniconda3/bin/teLeap: Warning!
Close contact of 0.988 angstroms between nonbonded atoms OD1 and H
------- .R<NASN 1>.A<OD1 11> and .R<LEU 2>.A<H 2>
/home/aiwt/miniconda3/bin/teLeap: Warning!
Close contact of 1.247 angstroms between nonbonded atoms HG and H
------- .R<LEU 2>.A<HG 9> and .R<TYR 3>.A<H 2>
...
注意: 这里的警告是分子结构中常见的非键原子距离过近警告. 使用vmd观察结构文件, 可以看到部分原子距离过近.
使用pymol查看多肽结构: pymol trp_initial.pdb, 从下图中可以看到明显的不合理键.
图. 1.5.1 肽链的初始结构#
1.5.3. 能量最小化#
肽链初始伸展构象可能存在原子近距离接触，需通过能量最小化去除应力，避免模拟初期能量异常。以下是能量最小化参数文件:
Minimization
&cntrl
imin=1, # 开启能量最小化
maxcyc=5000, # 总步数
ncyc=2500, # 前2500步采用最陡下降法，后2500步采用共轭梯度法
cut=10.0, # 非键相互作用截断距离（Å）
ntb=0, # 关闭周期性边界条件（使用于隐式溶剂模型）
igb=8, # 指定GB-Neck2隐式溶剂模型（igb=8对应GB-Neck2）
/
运行能量最小化: sander -O -i em.in -p trp.parm7 -c trp.crd -o em.out -r em.rst
提取结构文件: ambpdb -p trp.parm7 -c em.rst > trp_em.pdb
使用pymol程序查看多肽结构: pymol trp_em.pdb, 可以看到初始的不合理结构消失.
图. 1.5.2 经过能量最小化后的肽链结构#
1.5.4. 升温平衡#
将体系缓慢加热超过50 ps, 可分阶段进行, 分阶段对体系在每个温度下均衡加热可以降低系统崩溃的几率, 尤其是在处理你不熟悉的新的体系时, 这里分为两个阶段: 0 -> 300K; 300 -> 325K.
论文中提到有必要将系统加热到 325 K，以避免在动力学上陷入势能陷阱. 在高温运行时, 需要仔细考虑时间步长, 因为系统中的振荡会更加明显, 因此温度越高时间步长应该越小, 这里采用2fs足够.
(1). 第一阶段: 0 -> 300K
折叠模拟
&cntrl
imin=0, # 关闭能量最小化，开启MD
irest=0, #从初始结构开始新模拟（非重启）
ntx=1, # 仅读取初始坐标（不读取速度）
nstlim=50000, dt=0.0005,# 总步数5 ns，时间步长2 fs）
ntc=2, # 约束键长（SHAKE算法）
ntf=2, # 约束键长导数
ntt=1, # 采用 Berendsen 热浴进行温度耦合
tautp=1.0, # 温度耦合时间常数（1 ps）
tempi=0.0, temp0=300.0, # tempi 为初始温度，temp0 为目标温度
igb=8, # GB-Neck2隐式溶剂
ntb=0, # 无周期性边界（隐式溶剂无需盒子）
ntpr=50, # 每50步输出能量日志
ntwx=50, # 每50步保存轨迹
/
运行模拟: sander -O -i heat.in -p trp.parm7 -c em.rst -o heat.out -r heat.rst -x heat.nc
提取结构文件: ambpdb -p trp.parm7 -c heat.rst > trp_heat.pdb
(2). 第二阶段: 300 -> 325K
heating Stage 2
&cntrl
imin=0, irest=0, ntx=1,
nstlim=50000, dt=0.0005,
ntc=2, ntf=2,
ntt=1, tautp=1.0,
tempi=300.0, temp0=325.0,
ntpr=50, ntwx=50,
ntb=0, igb=8,
/
运行模拟: sander -O -i heat2.in -p trp.parm7 -c heat.rst -o heat.out -r heat2.rst -x heat2.nc
提取结构文件: ambpdb -p trp.parm7 -c heat2.rst > trp_heat2.pdb
图. 1.5.3 经过两次升温平衡后肽链结构变化#
可以看到肽链开始形成α螺旋, 距离稳定的折叠结构还需要继续模拟.
1.5.5. 成品 MD#
最后在 325 K 下进行长时间的成品模拟, 同样分阶段运行, 每次5ns, 使用下面的md.in参数文件, 每段结束后观察多肽的结构, 重复直至达到稳定结构.
md 1 0-5ns
&cntrl
imin=0, irest=1, ntx=5,
nstlim=2500000, dt=0.002,
ntc=2, ntf=2,
ntt=1, tautp=0.5,
tempi=325.0, temp0=325.0,
ntpr=500, ntwx=500,
ntb=0, igb=8,
/
运行命令:
sander -O -i md.in -p trp.parm7 -c heat2.rst -o md.out -r md.rst -x md.nc
ambpdb -p trp.parm7 -c md.rst > trp_md.pdb
sander -O -i md.in -p trp.parm7 -c md.rst -o md2.out -r md2.rst -x md2.nc
ambpdb -p trp.parm7 -c md2.rst > trp_md2.pdb
sander -O -i md.in -p trp.parm7 -c md2.rst -o md3.out -r md3.rst -x md3.nc
ambpdb -p trp.parm7 -c md3.rst > trp_md3.pdb
sander -O -i md.in -p trp.parm7 -c md3.rst -o md4.out -r md4.rst -x md4.nc
ambpdb -p trp.parm7 -c md4.rst > trp_md4.pdb
可以看到运行20ns后，结构达到稳定，与NMR结构类似。在pymol中使用super命令执行“分子叠加”可以对比两个结构的相似性。
图. 1.5.4 模拟过程中肽链的结构变化及NMR结构#
1.5.6. 结果分析#
进一步，可以对模拟体系的温度、势能、RMSD等进行分析，对模拟轨迹、氢键、簇分析等。