• 第1部分: 非标准残基的电荷生成
  • 1.1) 染料的电荷计算
    • i) 几何优化
    • ii) 计算静电势
    • iii) 将高斯ESP数据转换成RESP格式
    • iv) 为RESP生成输入文件
    • v) 运行染料多重构象的RESP拟合
  • 1.2) linker的电荷计算
    • i) 优化几何
    • ii) 计算静电电位
    • iii) 将高斯ESP拟合转换成RESP格式
    • iv) 为RESP生成输入文件
    • v) 运行RESP 拟合的阶段1
    • vi) 生成RESP的输入文件
    • vii) 运行RESP fit的第二阶段
• 第2部分
  • 2.1) linker和荧光素的连接过程
    • 2.1.1) 用Sirius改变荧光素和linker的坐标
    • 2.1.2) 删除剩余的优化封端基团并保存pdb文件
  • 2.2) 修改fam5.pdb文件
• 第3部分: 为FAM5构建一个Leap库文件
  • 3.1) 定义FAM5拓扑
  • 3.2.1) 定义原子类型和前面计算的原子电荷
  • 3.2.2) 定义一个不相似的氧原子的原子类型
  • 3.3) 保存新的库文件
  • 3.4) 创建一个frcmod文件.
• 第4部分 将FAM5分子连接到polyAT十聚体的5’末端
  • 4.1) 准备polyAT十聚体
  • 4.2) 使用Sirius将FAM和polyAT分子合并在一起
  • 4.3) 将组合的pdb加载到Leap并创建prmtop和inpcrd文件

• Bryan Leland, David Paul, Brent Krueger, Ross Walker(amber_tutorial_query@rosswalker.co.uk) 原始文档
• 2018-03-20 21:13:36 翻译: 陈亚玉(福州大学)

注意: 这些教程是为了提供如何使用AMBER软件套件来进行模拟的说明性例子, 这些模拟可以在合理的时间内在简单的工作站上运行, 但对特定的应用领域并不一定提供最佳的参数或方法选择.)

在本教程中, 我们将通过完成一个高级场景来准备一个体系并用sander来进行模拟, 这其中包括几个非标准残基, 在此用Antechamber 和GAFF力场可能不太合适(参见http://www.ambermd.org/tutorials/basic/tutorial4/index.htm使用Antechamber模拟的例子). 为了说明这个过程, 我们将选择一个实际的研究问题, 而不是一个简单的例子. 这当然有几个说明. 虽然这提供了一个更现实的教程, 但它当然存在任何研究项目中固有的问题. 也就是说, 我们通常有多种方法去完成某件事, 而我们并不一定需要事先知道哪一个是最好的选择. 正如你稍后看到的, 在推导某些RESP电荷的时候, 情况就是如此. 但是, 这种方法应该可以为你准备自己的体系时提供坚实的基础. 你只需知道, 当你的体系很复杂时, 没有教程可以一步一步指导你, 所以你应该把这个教程仅作为解决类似问题的基石. 而且大量的想法和文献搜寻对解决自己的体系也是必须的.

本教程主要是为荧光素染料结合到poly-AT DNA十倍体创建prmtop和inpcrd文件. 下面展示了这个拓扑的三维以及二维图.

图一

图二

DNA标注为蓝色, 荧光染料及其linker显示为红色. 这种设置通常用于激光探针实验. 在本教程中, 我们将通过一些步骤来创建prmtop和inpcrd文件. 这些步骤的大概顺序是:

1. 生成染料和连接基团的电荷.
2. 创建一个关于染料和linker的PDB文件, 将在步骤3使用.
3. 建立包含自定义参数的染料和linkerLeap库文件和frcmod文件.
4. 准备polyAT的5’末端用于连接染料和linker, 将染料linker连接到polyAT并保存prmtop和inpcrd文件.

注意: 本教程假定您可以访问某些不属于AMBER套件的程序. 这包括高斯(http://www.gaussian.com)和天狼星(http://www.ngbw.org/sirius/已失效). Sirius是免费的, 你可以下载并安装它, 但高斯是商业. 如果你没有Gaussian, 你仍然可以完成本教程, 因为提供了所有的输入和输出文件. 也可以使用免费提供的GAMESS包来计算RESP拟合的ESP, 而不是高斯. 但是, 这不在本教程中.

第1部分: 非标准残基的电荷生成

AMBER力场是围绕精确的成对电荷概念建立的. 电荷的基本单元是以比如以一个氨基酸为一个残基的形式. 在AMBER中, 每个残基被限定为具有一个单元的电荷(除了某些DNA和RNA末端残基). 这使得力场电荷模型可以转移, 因为单个残基可以用作构建模块来构建更大的蛋白质和核酸体系而不需要为每个模型重新调整电荷. 这种方法的关键在于电荷产生的方法是力场参数化的组成部分. 在AMBER中, 用于电荷产生的方法被称为约束静电势(RESP)方法. 有关RESP程序的细节和背景, 请参见J. Phys. Chem, 1993,97,10269-10280.

该程序的基础是确定旋转地简并的原子的多级方法, 如甲基和亚甲基氢具有相同的电荷. 这个过程的一个重要部分是这些成对电荷是从气相QM计算产生的, 但这些成对电荷被设计用来重现液相电荷的相互作用. 对于FF94, 96和99力场的变量, 是通过使用HF/6-31G* QM算法产生静电势来实现的. HF/6-31G*计算中的误差与气相和溶液中电荷分布之间的差异相似. FF03力场使用不同但功能类似的方法. 这里主要的要点是, 在产生电荷时, 不应该随意改变QM理论的层级.

电荷计算也是围绕一个具有完整电荷的残基的概念设计的. 理想状态下, 残基应该相当小, 以减少产生伪像的可能性, 因为QM计算或RESP拟合陷入高能量阈值. 因此在我们的例子中, 我们将染料和linker分成两个独立的残基, 这也有实际的原因. 例如, 我们可能想要使用许多带有相同linker的不同染料.

染料

链接

图1.1

AMBER中有多种电荷拟合可供选择. 这包括最初的RESP执行, Antechamber和R.E.D. 对于独立的配体, 通常使用Antechamber(参见http://www.ambermd.org/tutorials/basic/tutorial4/index.htm, 一个用Antechamber生成坐标文件和力场文件的例子), 但它不是用来处理共价结合的残基, 当然你可以调整某些参数去做, 但它对关于旋转的简并化原子仍存在问题. 目前对于大多数复杂的系统来说, 最好的选择就是使用R.E.D. 因为它提供了更多的重复性, 而且是用来自动产生多方向拟合等, 也使其他人更容易重现电荷计算.

由于我们的意图是重现这个特定的研究项目所使用的内容, 并强调了RESP拟合所涉及的实际程序, 所以我们将使用RESP手动拟合电荷. 值得注意的是, 这里所使用的程序和所选的选项是征对于这个研究项目的. 这些程序并不是逐步指导在AMBER力场产生”标准”电荷. 这仅仅是为了说明这涉及的典型步骤. 欲了解更多有关原始力场如何拟合和程序如何应用于未来的力场, 请参考文献.

1.1) 染料的电荷计算

电荷拟合由3部分组成. 首先我们需要优化几何结构, 然后我们需要产生一个静电势, 最后我们需要运行RESP程序. 然而, 在我们做之前, 我们通常需要处理染料与linker连接的键. 在开始几何优化之前, 这需要封端. 理论上可以使用任何形式的封顶, 但是, 重要的是要记住封端原子不会出现在实际的模拟中, 而非封端原子上的电荷必须加到正确的整数电荷上. 对于蛋白质, 作为力场设计的一部分而开发的封端程序是使用ACE或NME残基:

ACE

NME

图1.2 改编自Cornell等人JACS 117, 1195, 5179

上图显示了在AMBER FF9X力场中用到的两个封端基团. 你应该注意到, 这两个封端基团的电荷总和为零. 我们将用NME将染料的共价键断开, 而后执行RESP程序时, 我们将限制NME原子使其与FF9X力场中的电荷相同.

i) 几何优化

第一步是优化已经封端的染料. 这需要在相当高的层级里用QM方法来完成. 最初的AMBER FF94中用MP2/6-31G* 来优化几何构型. 在这里为了速度, 我们将使用B3LYP/6-31G*, 它应该与MP2的结果相似. 你可以在您选择的QM软件中进行优化. 在这种情况下, 我们使用高斯. 高斯输入文件的创建超出了本教程的范围, 但下面提供了参考. 结果显示, 染料在不同的构象中有两个极小值, 且能量相当. 我们称这两个构象为floF和floB, 两者如下所示:

floB.pdb

floF.pdb

图1.3 注意Amide组的方向.

为了处理这种情况, 我们将做一个多构象拟合. 第一步就是优化两个结构.

这里是高斯输入文件: floB_opt.gin, floF_opt.gin

这里是输出文件: floB_opt.gout, floF_opt.gout

我们还将生成一个优化后的floB结构的pdb以供在leap中构建单元. 您可以使用多种方法从高斯生成pdb. 带有Gaussian的Newzmat可以将一个检查点文件转换成一个pdb文件, 但它非常麻烦, 而且经常崩溃. 另一个选择是使用Gaussview. WebMO(http://www.webmo.net/)也可以用来设置和运行高斯计算, 你也可以使用Babel. 本教程在这里提供了pdb: floB_opt.pdb.

ii) 计算静电势

下一个步是计算这两个结构优化过的ESP以便RESP程序读取. 这的关键在于我们使用HF/6-31G*作为理论水平由于如前所述随机消除错误. 这个工作的高斯路线需要读:

#P HF/6-31G* SCF=Tight Pop=MK IOp(6/33=2, 6/41=10, 6/42=17)

• IOp 6/33=2: 使高斯输出势能点和电位 - 不要去改动.
• IOp 6/41=10: 指定每个原子使用10个同心层的点.
• IOp 6/42=17: 指定每层中的密度点. 17 赋予每个原子2500个点. 对于较大的分子, 可能需要稍微减少点数. RESP可以用最多的ESP点数是99,999. 17的值可以减小到任何整数值.

输入文件是从优化步骤中生成的最终结构. 注意这里使用6/41=10和6/42=17是该项目为增加高斯产生的格点数. 对于大多数”标准”RESP拟合只需要IOp(6/33=2).

注意我从优化中得到检查点文件并将它们复制到floB_hf.chk和floF_hf.chk中, 这样高斯输入文件就可以只读geom = check而不覆盖原来的检查点文件.

这里是输入文件: floB_hf.gin, floF_hf.gin

这里是输出文件: floB_hf.gout, floF_hf.gout

iii) 将高斯ESP数据转换成RESP格式

RESP只读取特定文件格式的ESP数据. 要从高斯输出格式转换为RESP输入格式, 我们将使用以下脚本:

1
2
3
4
5
6
7

#!/bin/csh
xlf /usr/local/apps/amber9/src/resp/readit.f
grep "Atomic Center " $1 > a
grep "ESP Fit" $1 > b
grep "Fit    " $1 > c
./a.out
rm -f a b c a.out readit.o

注意: 您需要根据你的电脑修改这个脚本. 要使用正确的编译器(例如g77)和readit.f的正确路径. 如果您在Cygwin中运行此文件, 那文件名称要为a.exe或a.out.exe, 所以请适当的更新上述脚本. 另外请注意在AMBER 10中espgen替代了readit.f, 但是这个教程还没有更新, 所以如果你使用的是AMBER10, 你就需要用这个链接下载readit.f文件.

为了运行这个脚本, 我们首先需要确定有多少个原子, 以及有多少个ESP中心. 这是ESP数据生成之前的高斯输出文件. 例如:

Atomic Center 1 is at 7.127640 .120744 -.032564
Atomic Center 2 is at 7.345002-0.099693 -1.079630
...
Atomic Center 42 is at -4.823702 2.342669 -.413287
ESP Fit Center 43 is at 7.127640 .120744 2.067436
...
ESP Fit Center**** is at -6.585767 1.848962 -2.848473
ESP Fit Center**** is at -6.502128 1.613625 -2.848473

所以原子计数很简单. 这里是42. ESP拟合中心有点复杂因为它超过9999和高斯在这就只输出*s. 所以你知道行号之间的区别(例如使用一个程序, 如nl ‘ nl floB_hf.gout> floB_hf.gout_line_numbers’). 在floB例子中, 我们有94,681个ESP拟合中心, 而 floF有94,998个.

我们现在可以运行esp.sh:

>./esp.sh floB_hf.gout
enter natom,nesp:

所以我们输入42, 94681.

这将会输出一个esp.dat文件, 我们将其重命名floB_esp.dat.

然后重复这个做floF, 你最终应该有两个dat文件: floB_esp.dat, floF_esp.dat

最后我们需要把这两个dat文件合并起来, 因为RESP需要所有的数据在一个文件中进行多重构象拟合floBF_esp.dat:

> cat floB_esp.dat floF_esp.dat> floBF_esp.dat

iv) 为RESP生成输入文件

为了做RESP拟合, 我们将遵循Cornell等人在最先RESP文献中使用的相同方法.

Cornell, W.D., Cieplak, P., Bayly, C.I., Kollman, P.A., JACS, 1993, 115, 9620-9631. (pdf)

通常执行RESP拟合需要两个步骤. 首先所有的原子都允许被变化, 接着就是所有的简并氢原子被定为具有相同的电荷和自拟合. 然而染料没有甲基或亚甲基(除了已经NME甲基基团), 因此它们不需要执行第二步.

我们还需要记住, 我们要做一个多重构象的RESP拟合, 因为我们有两个相似的状态.

这是第一步的输入文件: floBF-resp.in

你需要认真看这个文件以确保你理解它. 第一步包括:

floBF-resp run #1
 &cntrl
 nmol=2,
 ihfree=1,
 qwt=0.0005,
 iqopt=2,
 /

nmol = 2: 告诉RESP, 这是一个适合超过两个”分子”(染料的两个构象)的多重构象拟合. ihfree = 1: 告诉RESP我们只想要对重原子的弱约束. qwt = 0.0005: 约束强度是0.0005 A.U. iqopt = 2: 告诉RESP以拉格朗日约束的形式从qin文件读取初始电荷. 因此这是我们可用来修改NME封端电荷.

下一部分包含两个分子:

	1.0
floB
   -2 42
	6 -1
	1 -1
	1 -1
	1 -1
	7 -1
	1 -1
	6  0
	8  0
	6  0
	6  0
	1  0
	6  0
	6  0
	8  0
	8  0
	6  0
...

第一个值(1.0)是多分子拟合的加权因子. 在这种情况下, 我们将它们设置为1.0, 因为我们希望它们的权重相等. 下一行(floB)只是一个标识符注释. 然后我们有一条线, 指定总电荷和原子数(-2 42). 以下每行指定与QM计算中相同顺序列出的单个原子. 第一列是原子序数, 第二列是N.

N可以是-1,0或正整数. 如果N = -1, 那么它告诉RESP这个原子的电荷应该从qin文件中读取并且被约束到这个值. 0表示可以在拟合过程中自由更改, 正整数表示允许你指定两个原子具有相同的电荷. 这用于拟合的第二步, 以约束简化原子具有相同的电荷. 例如, 如果我们希望原子11与原子10具有相同的电荷, 但要使它们成对地变化, 那么原子10将为0, 原子11为10.

最后一部分是做多重构象拟合的.

	2
	1 1 2 1
	2
	1 2 2 2
	2
	1 3 2 3
	2
	1 4 2 4
	2
	1 5 2 5
	2
	1 6 2 6
	2
	1 7 2 7
	2
...

在这里你需要为每个原子设置2行. 第一行只是指定分子的数量, 所以在这种情况下它都是2.第二行的格式为:

分子 原子 分子 原子

并用于指定分子之间的等价性. 所以在这里第一行读: 1 1 2 1.这意味着分子1的原子1等同于分子2的原子1(‘ 我们知道在这个例子中这看起来很愚蠢, 但是在某些情况下它可能是有意义的. ‘)

接下来我们需要一个qin文件来为封端提供电荷.

floBF-resp.qin

 -0.149000 0.097600 0.097600 0.097600 -0.415700 0.271900 0.000000 0.000000
  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
  0.000000 0.000000 -0.149000 0.097600 0.097600 0.097600 -0.415700 0.271900
  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
  0.000000 0.000000 0.000000 0.000000

每行有8个指定每个分子中的每个原子的电荷的条目, 其与QM计算中相同的顺序. 在这种情况下, 我们设置所有原子, 除了NME原子的初始电荷我0. 我们设定的NME原子具有来自FF99SB力场(与FF99 / 94力场相同)的电荷. 注意那些不为零的条目在RESP输入文件中必须等于-1.

v) 运行染料多重构象的RESP拟合

> $ AMBERHOME / exe / resp -O -i floBF-resp.in -o floBF-resp.out -p floBF-resp.pch -t floBF-resp.chg -q floBF-resp.qin -e floBF_esp.dat

输入文件: floBF-resp.in, floBF-resp.qin, floBF_esp.dat

输出文件: floBF-resp.out, floBF-resp.pch, floBF-resp.chg

输出文件的检查将显示每个原子的电荷如下:

图1.4

如上所述, 除了始终保持固定的封端以外, 我们不需要进行第二步的配合, 因为这里没有甲基或亚甲基. 我们也不需要担心COO组中的简并性, 因为这里的文献搜索显示这里的旋转屏障在6kCal/mol左右, 这意味着它们不一定是简化的.

你应该记下这些电荷, 因为我们以后需要这些电荷.

1.2) linker的电荷计算

我们现在重复图1.1所示的linker的过程. 然而如何限制DNA末端的linker有点棘手. 染料末端容易, 我们只是使用ACE, 但DNA很难. 还有一个问题, 即封端被移除后的linker必须有-0.3079的电荷. 其原因是通常DNA有5’和3’端. 在AMBER FF9 9力场中, 标准的”内部”dna(和rna)残基有-1的整数电荷. 然而, 末端碱基成对参数化, 当加在一起时给出-1的整数电荷. 5’端残基的电荷为-0.3079, 而3’端的电荷为-0.6921. 在连接linker与DNA时, 我们将用”内部”的残基取代原来的5’残基. 因此, 现在linker是有效的5’残基, 所以负责弥补剩余的部分电荷. 因此, linker上的DNA封端基团需要被约束为-0.6921. AMBER力场中没有残基有类似的性质, 因此我们需要”发明”我们自己的封端基团. 这里有很多不同的选择. 理想情况下, 人们会创建几个不同的封端基团, 并尝试所有这些来看看最终的电荷是如何变化的. 在这里, 我们只提出一个可能的, 但绝不是完美的解决方案. 图1.5(左图)显示了我们的linker在一端被ACE封端, 另一端连接到磷酸盐和糖, 因为它将在我们的最终DNA系统中. 右边是我们将在RESP拟合中使用的近似值. 因此, linker现在是有效的5’端残基, 所以负责弥补剩余的分数电荷. 因此, linker上的DNA封端基团需要被约束为-0.6921. AMBER力场中没有残基, 因此我们需要”发明”我们自己的封端基团. 这里有很多不同的选择. 理想情况下, 人们会创建几个不同的封端基团, 并尝试所有的封端基团来看看最终的电荷是如何变化的. 在这里, 我们只提出一个可能的, 但绝不是完美的解决方案. 图1.5(左图)显示了我们的linker一端被ACE封端, 另一端连接到最终DNA系统中的磷酸盐和糖. 右图是我们将在RESP拟合中使用的近似结构.

由于linker与DNA残基上的磷酸键结合, 因此使用磷酸盐作为封端看似挺合理的. 这个想法与教程中描述的在AMBER力场中计算DNA的电荷伴着计算磷酸盐的甲基端一样.

Cieplak, P., Cornell, W.D., Bayly, C., Kollman, P.A., J.Comp.Chem., 1995, 16(11), 1357-1377 (pdf)

图1.5

我们需要限制封端的原子电荷总和为-0.6921, 但是这样做的方式是尽可能保留能代表DNA残基的电荷场. 我们首先进行一个RESP拟合, 就像我们对染料-NME结构所做的那样, ACE的电荷设定在FF99SB力场中(见图1.2). 对磷和三个封闭基团的氧原子(图1.6中以蓝色突出显示)的电荷与FF99SB DNA电荷的电荷相似(见图1.7). 封端基团(PO 2 -O-CH 3)被约束为总电荷等于-0.6921. 由于在FF99SB力场中没有等价基团可以修正它们, 因此允许甲基自由拟合(在基团约束内). 由于ACE基团是中性的, 并且是整个构造的总电荷是-1, 我们的linker的电荷为-0.3079. 为了进行比较, RESP拟合没有限制地运行, 所得到的电荷是相似的.

图1.6

图1.7 改编自Cornell等人 JACS 117,1195,5179

如前所述, 我们首先优化几何形状, 然后在进行RESP拟合之前计算静电势. 这些文件在下面提供.

i) 优化几何

高斯输入文件: ln5_opt.gin

输出文件: ln5_opt.gout

PDB文件: ln5_opt.pdb

在运行HF ESP计算之前, 将优化运行中生成的ln5_opt.chk复制到ln5_hf.chk, 以便可以使用geom = check.

ii) 计算静电电位

高斯输入文件: ln5_hf.gin

输出文件: ln5_hf.gout

iii) 将高斯ESP拟合转换成RESP格式

1
2
3
4
5
6
7

#!/bin/csh
xlf /usr/local/apps/amber9/src/resp/readit.f
grep "Atomic Center " $1 > a
grep "ESP Fit" $1 > b
grep "Fit    " $1 > c
./a.out
rm -f a b c a.out readit.o

再一次, 请记住根据需要编辑这个脚本.

确定有多少个原子, 以及有多少个ESP拟合中心: 原子= 35, ESP拟合中心= 79216.

运行esp.sh:

>./esp.sh ln5_hf.gout
enter natom,nesp:

输入35, 79216.

这将会输出一个esp.dat文件, 我们将重命名ln5_esp.dat.

由于我们在这里只做一个单一的结构拟合, 所以在这里不需要合并任何文件.

iv) 为RESP生成输入文件

接下来, 我们生成RESP拟合所需的输入文件. 这些与我们用于上述染料的那些类似, 但是包含一些细微的变化, 这些变化仅涉及单个构象拟合和使用基团电荷约束.

这是第一阶段的输入文件: linker_resp.in

ln5 RESP run #1
 &cntrl
 ihfree=1,
 qwt=0.0005,
 iqopt=2
 /

这一次我们没有指定nmol, 因为我们将其默认为1, 因为我们不再进行多重构象拟合.

下一部分只包含一个分子, 如上所述, 我们将ACE和选定的磷酸盐封端原子固定(-1)并从qin文件读取:

	1.0
Link
   -1 35
	6 -1
	1 -1
	1 -1
	1 -1
	6 -1
	8 -1
	7 0
	1 0
	6 0
	1 0
	1 0
	6 0
	1 0
	1 0
	6 0
	1 0
	1 0
	6 0
	1 0
	1 0
	6 0
	1 0
	1 0
	6 0
	1 0
	1 0
	8 0
   15 -1
	8 -1
	8 -1
	8 -1
	6 0
	1 0
	1 0
	1 0

我们不需要从前一个文件的最后一部分, 因为我们没有做多重构象拟合. 但是, 我们确实需要额外的部分来定义上面讨论的基团约束.

8 -0.69210
1 28 1 29 1 30 1 31 1 32 1 33 1 34 1 35

这里, 8是约束中的原子数. -0.69210是这些原子的电荷总和. 第二行指定基团中涉及的每个原子的分子数和原子数. 注意RESP预计每行不超过16个值, 所以如果你有超过8个原子, 这个基团描述将跨越多行.

接下来我们需要一个qin文件为封顶基团提供电荷.

linker - resp1.qin

 -0.366200 0.112300 0.112300 0.112300 0.597200 -0.567900 0.000000 0.000000
  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
  0.000000 0.000000 0.000000 1.165900 -0.776100 -0.776100 -0.495400 0.000000
  0.000000 0.000000 0.000000

v) 运行RESP 拟合的阶段1

>$AMBERHOME/exe/resp -O -i linker_resp.in -o linker_resp.out -p linker_resp.pch -t linker_resp.chg -q linker-resp1.qin -e ln5_esp.dat

输入文件: linker_resp.in, linker-resp1.qin, ln5_esp.dat

输出文件: linker_resp.out, linker_resp.pch, linker_resp.chg

首先, 我们可以看一下输出文件(linker_resp.out), 以检查阶段1优化的电荷是否合理, 受约束的电荷是否保持不变, 是否遵守基团约束.

vi) 生成RESP的输入文件

我们现在需要准备RESP fit的第二阶段(linker_resp2.in). 在这里, 我们将保持除甲基和亚甲基以外的所有原子电荷不变. linker几乎是整个分子, 除了NH和O之外的所有其他元素. 除了指定附加原子被固定外, 我们也对输入文件做了一些小的修改. 第一部分现在包含:

ln5 RESP run #2
 &cntrl
 ihfree=1,
 qwt=0.001,
 iqopt=2
 /

除了我们将权重(qwt)更改为0.001之外, 这与第一次运行相同.

文件的第二部分比第一阶段稍微复杂一点. 在这里, 我们把所有不是甲基或亚甲基里氢/碳的基团都固定(-1)并从qin文件中读取. 甲基/亚甲基碳, 每个基团的第一个氢我们设置为零, 以便他们可以自由改变. 我们设定每个基团中加入的氢的电荷都受到第一个氢的约束. 这是通过为在第二列中每个额外的氢中指定第一个氢的原子数来完成的. 为了保持一致, 我们还将重新拟合我们的DNA上的封端基团, 但保持基团的电荷为-0.69210.

	1.0
Link
   -1 35
	6 -1
	1 -1
	1 -1
	1 -1
	6 -1
	8 -1
	7 -1
	1 -1
	6 0
	1 0
	1 10
	6 0
	1 0
	1 13
	6 0
	1 0
	1 16
	6 0
	1 0
	1 19
	6 0
	1 0
	1 22
	6 0
	1 0
	1 25
	8 -1
   15 -1
	8 -1
	8 -1
	8 -1
	6 0
	1 0
	1 33
	1 33
	8 -0.69210
   1 28 1 29 1 30 1 31 1 32 1 33 1 34 1 35

vii) 运行RESP fit的第二阶段

现在我们可以运行RESP的第二阶段了. 注意我们使用在阶段1中的输出文件linker-resp.chg作为阶段2的qin(-q)文件.

>$AMBERHOME/exe/resp -O -i linker_resp2.in -o linker_resp2.out -p linker_resp2.pch –t linker_resp2.chg -q linker_resp.chg -e ln5_esp.dat

输入文件: linker_resp2.in, linker_resp.chg, ln5_esp.dat

输出文件: linker_resp2.out, linker_resp2.pch, linker_resp2.chg

现在我们可以检查从RESP拟合的第二阶段生成的输出文件, 并因此标记所有原子的所有电荷. 这在下图中总结, 尽管最终我们将得到所有关于linker原子的电荷:

图1.8

我们现在准备进行到第二部分, 我们将创建linker和染料组合的pdb文件.

第2部分

2.1) linker和荧光素的连接过程

linker作为荧光素的锚点, 使我们能够将荧光染料附着到polyAT十聚体的A5’末端. 因为这样做, 我们能够模拟许多激光实验中存在的探针连接. 我们将使用 Sirius可视化程序(关于Sirius的教程就在这里), 尽管原则上你可以使用你另一个程序. 你可以在xleap中做, 但这有点困难. 这一步是必需的, 因为没有这个染料, linker和DNA系统的晶体结构. 因此, 我们必须自己建立一个合适结构的PDB文件, 我们最终会加载这个文件当在在生成prmtop和inpcrd文件时. 如果要模拟的系统的晶体/核磁共振结构是可用的, 那么根据需要编辑这个PDB文件就是一件简单的事情, 本节讨论的大部分方法就不再需要了.

由于我们在上一节中分别对荧光素染料和连接分子进行了优化, 因此我们必须将这两个分子键合在一起, 正如引言所描述的(如图2所示). 幸运的是, 这有几种工具可以让我们操作分子并进行必要的改变. Sirius可视化程序使我们能够将荧光素的NH基团(NME封端)叠加在存在于linker中的NH基团上, 以及两个分子中存在的羰基上. 然后删除这些”封端”, 并将荧光linker坐标保存在相同的参考框架中, 从而可以将这些分子在xleap里连接在一起.

考虑到这些, 执行Sirius程序, 并从先前电荷拟合(floF_opt.pdb, ln5_opt.pdb)中加载这两个优化结构. Select, Tools > Structure Browser以显示这两个文件的内容. 这将派上用场, 因为它使得整个分子和原子片段的选择变得容易.

图2.1

2.1.1) 用Sirius改变荧光素和linker的坐标

为了选择我们要叠加的原子, 我们需要将linker从荧光素分子上移开. 首先, 在Structure Browser中”加亮”linker, 并选择位于工具栏中的”Set separate motion”按钮.

鼠标左键可旋转选中的原子, 中间的按钮放大, 右键平移原子. 有时, 旋转两个分子以获得更好的视角是有作用的. 这可以通过选择位于工具栏中的”Set common motion”按钮来完成. 注意: Mac OSX用户可能无法查看按钮工具提示, 因为Sirius仍处于测试阶段. 这个错误应该在更高版本的软件中修复. 有关Sirius按钮及其功能的完整列表, 请选择帮助> Sirius帮助…, 然后使用索引选项卡搜索”工具栏”.

图2.2

现在我们可以区分这两个分子, 我们可以重叠封端基团并删除额外的原子. 选择, 结构>结构叠加>原子锚定叠加. 应该出现下面的菜单(图2.3).

你可以使用鼠标来选择重叠过程中使用的原子. 这些原子以两个一组(一个是对准结构, 另一个是参考结构)选择. 基于本教程的目的, 哪个结构是”对齐的”, 哪个是”参考”并不重要. 因此, 我们已经任意地将linker作为”对齐”结构. 注意在每对中选择的第一个原子必须来自对齐的结构, 而第二个来自参考结构. 从linker中选择第一个原子N, 从染料中选择第二个原子N. 重复此过程 选择羰基原子C和O, 然后按OK. 由此产生的结构应该类似于图2.4的结构.

图2.3

图2.4

2.1.2) 删除剩余的优化封端基团并保存pdb文件

我们将称为组合的荧光素/linker分子为FAM, FAM几乎是完整. 下一步将按本教程前面讨论的封端基团删除. ACE和NME封端仅用于电荷计算, 并且在最终的FAM分子中不存在额外的原子. 荧光素的NME封端含有以下原子: CH3NH(图1.2). 使用Structure Browser突出显示这些原子, 并选择Structure> Edit Structure Objects. 一个GUI会出现; 确保选择”delete a set of atoms”按钮, 以及”currently selected atoms”按钮, 然后按OK. 对于linker上的ACE封端, 应重复上述过程. ACE封端上的原子是CH3CO(图1.2). 最后, 磷酸封端基团也需要从linker中去除. 选择POOO-CH3原子并删除. 所有额外的原子被删除后, 结构应该类似于图2.5的结构.

图2.5

我们现在已经用Sirius成功地修改了我们的分子, 并且准备好继续下一步的教程. 但是, 为了使进程顺利进行, 我们需要将我们的结构合并到一个pdb中, 而后编辑这个文件, 这样就可以正确被读入到xleap程序中. Sirius当前版本中的命名传统不允许使用3位以上的原子名称, 这对于使用荧光素这样的大分子的研究而言是个问题. 其次, 合并结构有两个不同的残基名称. 我们将创建一个库文件, 告诉xleap如何识别荧光素/linker结构, 因此我们的pdb中的每个原子必须在相同的残基中, 并且必须是相同的名称. (注意原则上我们也可以按照这两个残基设置两个库文件. )第三, 在新的pdb文件中的连接信息没有用, 所以需要删掉. 注意, 其中几个问题将在Sirius的更新版本中得到修复, 因为它仍处于测试阶段. 在此之前, 我们将手动解决这些问题.

所以, 我们可以继续合并我们的分子到一个pdb中. Sirius可以为我们自动完成这一步, 并减少需要完成的文本编辑工作. 要完成此任务, 请选择 Structure > Edit selected objects . 确保你选择”merge two structures together”按钮. 接下来, 我们将使用提供的下拉菜单来选择荧光素作为第一个结构, linker作为第二个结构. 设置这些变量后, 我们将命名新的结构fam5, 然后单击确定. 合并结构后, 我们可以使用File> Save结构来保存新的pdb文件. fam5.pdb

2.2) 修改fam5.pdb文件

在本教程的这一部分, 我们将修改新创建的fam5.pdb文件, 使其与xleap兼容. 这可以在您选择的文本编辑器中完成. 第一件事是编辑原子名称. 每个原子必须拥有一个自己的名字. 因为这是一个高级教程, 所以我们不会详细介绍pdb文件格式的标准原子命名方案. 接下来, 我们将为整个分子指定FAM的残基名称, 并将linker的链值由2改为1, 以使其被认为是第一个残基的一部分. 对pdb的最后一项主要改动是删除Sirius提供的连接信息. 我们从pdb文件中得到的数据是我们的荧光素和linker原子的坐标. 修改后的pdb应如下面的文件所示: fam5_edit.pdb. 我们现在准备将编辑后的pdb加载到xleap中, 并开始构建我们的库文件.

第3部分: 为FAM5构建一个Leap库文件

本教程的第三部分将介绍如何使用第一部分中推导出的电荷为我们的非标准残基(FAM5)构建一个leap库文件. 本教程中使用的FAM5残基在两个分块中进行了优化, 然而, 因为我们使用Sirius可视化程序”combined”这些分块, 我们只需要为整个FAM5分子构建一个库文件. 这个库文件将包含在我们的系统上运行MD所需的所有参数和电荷. 以下步骤的简要概述包括将FAM5加载到xleap程序中, 定义结构的拓扑和原子类型, 以及将以前产生的电荷输入到库文件中.

3.1) 定义FAM5拓扑

因为我们在pdb中删除了连接信息, 所以我们需要手动告诉xleap什么原子被连接在一起. 幸运的是, 使用xleap的编辑GUI可以轻松实现这一点, 这使我们可以拉这键合的结构. 我们将需要指定我们pdb文件的目录并执行xleap程序.

$AMBERHOME/exe/xleap &

下面的GUI会出现:

图3.1

将fam5_edit.pdb文件加载到xleap:

FAM5 = loadpdb"fam5_edit.pdb"

屏幕上会出现以下输出:

Loading PDB file: ./fam5_edit.pdb
Unknown residue: FAM5   number: 0   type: Terminal/last
..relaxing end constraints to try for a dbase match
  -no luck
Creating new UNIT for residue: FAM5 sequence: 1
Created a new atom named: C within residue: .R<FAM5 1>
Created a new atom named: O within residue: .R<FAM5 1>
Created a new atom named: C1 within residue: .R<FAM5 1>
Created a new atom named: C2 within residue: .R<FAM5 1>
Created a new atom named: H2 within residue: .R<FAM5 1>
Created a new atom named: C3 within residue: .R<FAM5 1>
Created a new atom named: C4 within residue: .R<FAM5 1>
Created a new atom named: O41 within residue: .R<FAM5 1>
Created a new atom named: O42 within residue: .R<FAM5 1>
Created a new atom named: C5 within residue: .R<FAM5 1>
Created a new atom named: C6 within residue: .R<FAM5 1>
Created a new atom named: H6 within residue: .R<FAM5 1>
Created a new atom named: C7 within residue: .R<FAM5 1>
Created a new atom named: H7 within residue: .R<FAM5 1>
Created a new atom named: C8 within residue: .R<FAM5 1>
Created a new atom named: C9 within residue: .R<FAM5 1>
Created a new atom named: C10 within residue: .R<FAM5 1>
Created a new atom named: H10 within residue: .R<FAM5 1>
Created a new atom named: C11 within residue: .R<FAM5 1>
Created a new atom named: H11 within residue: .R<FAM5 1>
Created a new atom named: C12 within residue: .R<FAM5 1>
Created a new atom named: O12 within residue: .R<FAM5 1>
Created a new atom named: C13 within residue: .R<FAM5 1>
Created a new atom named: H13 within residue: .R<FAM5 1>
Created a new atom named: C14 within residue: .R<FAM5 1>
Created a new atom named: O15 within residue: .R<FAM5 1>
Created a new atom named: C16 within residue: .R<FAM5 1>
Created a new atom named: C17 within residue: .R<FAM5 1>
Created a new atom named: C18 within residue: .R<FAM5 1>
Created a new atom named: H18 within residue: .R<FAM5 1>
Created a new atom named: C19 within residue: .R<FAM5 1>
Created a new atom named: H19 within residue: .R<FAM5 1>
Created a new atom named: C20 within residue: .R<FAM5 1>
Created a new atom named: O20 within residue: .R<FAM5 1>
Created a new atom named: C21 within residue: .R<FAM5 1>
Created a new atom named: H21 within residue: .R<FAM5 1>
Created a new atom named: N within residue: .R<FAM5 1>
Created a new atom named: H within residue: .R<FAM5 1>
Created a new atom named: C22 within residue: .R<FAM5 1>
Created a new atom named: H221 within residue: .R<FAM5 1>
Created a new atom named: H222 within residue: .R<FAM5 1>
Created a new atom named: C23 within residue: .R<FAM5 1>
Created a new atom named: H231 within residue: .R<FAM5 1>
Created a new atom named: H232 within residue: .R<FAM5 1>
Created a new atom named: C24 within residue: .R<FAM5 1>
Created a new atom named: H241 within residue: .R<FAM5 1>
Created a new atom named: H242 within residue: .R<FAM5 1>
Created a new atom named: C25 within residue: .R<FAM5 1>
Created a new atom named: H251 within residue: .R<FAM5 1>
Created a new atom named: H252 within residue: .R<FAM5 1>
Created a new atom named: C26 within residue: .R<FAM5 1>
Created a new atom named: H261 within residue: .R<FAM5 1>
Created a new atom named: H262 within residue: .R<FAM5 1>
Created a new atom named: C27 within residue: .R<FAM5 1>
Created a new atom named: H271 within residue: .R<FAM5 1>
Created a new atom named: H272 within residue: .R<FAM5 1>
Created a new atom named: O3* within residue: .R<FAM5 1>
  total atoms in file: 57

发生了什么

Xleap在它的数据库中不能识别FAM5残基, 因此它开始创建一个”UNIT”或新的残基来定义这个未知的分子. 作为这个过程的一部分, xleap将所有的原子添加到UNIT中. 因为我们正在处理一个非标准的残基, 这正是我们想要的. 一旦我们保存我们的库文件, 由xleap创建的UNIT将使程序识别FAM5. 仔细检查xleap”创建”了正确的原子是很重要的. 如果有差异, 你应该检查你加载的pdb文件. 最容易发现的差异是文件中的原子总数, FAM5的值应该是57.

我们现在已经把荧光素分子加载到了Xleap中. 但是, 如前所述, fam5_edit.pdb文件只含有原子坐标的信息. 因此, xleap没有关于荧光素参数和键连的信息, 我们需要使用图形界面的xleap手动输入这些信息.

使用xleap执行编辑界面

edit FAM5

在开始之前, 您应该知道如何操作编辑界面. 鼠标左键控制原子选择, 鼠标中键旋转分子, 鼠标右键平移. 如果MIDDLE按钮不起作用, 可以通过按住Ctrl键和LEFT鼠标按钮来旋转它. 若要使图像变大或变小可同时按住MIDDLE和RIGHT按钮(或Ctrl键和鼠标右键)和向上移动鼠标放大或向下以缩小. 为了帮助你识别: 碳原子是绿色的, 氧红, 氢白和氮蓝. 当前选中的原子以紫色显示(您可以通过在空白地方单击来取消选择原子, 或者在弹性连接/选择时按住shift键). 您可以使用擦除工具撤消绘图过程中意外创建的键/原子. 此外, 您可以通过选择显示 - >名称来查看原子名称. 一旦你进入”绘制”模式, 你可以将原子键合在一起. 要连接原子, 请左键单击第一个原子, 然后不松开, 将光标拖到第二个原子上. 当逐渐靠近的时候, 键应该”吸附”到预期的位置. 松开鼠标时完成键的连接. 小心…注意到使用xleap时, 不要用X按钮关闭窗口. 此操作将终止会话, 任何未保存的数据将丢失. 这将是特别麻烦的当我们需要输入的大量的数据时. 我们应该是使用File - >close方法. 对于主编辑窗口, 可以使用Unit - > Close来完成.

图3.2

图3.3

3.2.1) 定义原子类型和前面计算的原子电荷

现在我们已经定义了FAM拓扑, 现在到处理库文件的时候了. 如前所述, 我们需要定义分子中每个原子的原子类型. 在力场文件和AMBER 数据路径里的文件中描述了ff99SB 力场中支持的原子类型列表. 我们将通过类比来定义我们的原子类型, 比较教程中的原子与我们实际系统中的原子. 这种类比将考虑到每个原子的键合和二面角相互作用. 在$ AMBERHOME / dat / leap / parm中的parm99.dat文件在描述每个原子类型的环境方面做得很好, 可以让我们做出准确的比较. 当我们处理这个过程时, 我们将会遇到一个没有被任何parm99原子类型所定义的氧分子. 因此, 我们将创建一个新的原子类型来描述这种氧的键/二面角相互作用, 并在库中定义它.

要定义这些参数, 我们将使用xleap中给出的原子描述表. 双击编辑窗口的黑色区域, 选择整个FAM分子. 一旦选定了所有的原子, 选择Edit> Edit selected atoms. 显示窗口就会出现这个表. 注意: Mac用户如果在查看整个”编辑选定原子”窗口时遇到问题, 应从X11工具栏中选择”窗口”>”缩放”.

图3.4

我们将在这个窗口中输入我们所有的原子类型和电荷. 由于我们系统中有大量的原子, 所以我们在教程中不会把每一个原子拿来做比较而是聚焦在几个原子; 给你选择自己分析剩下的原子. 在本教程中, 我们将向您提供包含所有的原子类型和电荷的一个文件, 可供你输入到”编辑选定的原子”窗口中.

原子比较 #1

在我们的第一个原子比较中, 我们将看一下C10, 它是fam5_edit.pdb文件中的17号原子. 该原子是氧杂蒽环的一部分, 类似于苯环中存在的碳原子, 因此我们指定了适合于芳香碳的原子类型CA.

原子比较 #2

在我们的第二个比较中, 我们将看到fam5_edit.pdb文件中第18原子H10. 该原子是键合到芳香碳上的氢原子. 因此, 我们将找一个适合于这种类型的氢原子类型, HA.

原子比较 #3

对于第三个例子, 考虑原子O41和O42. 这些是在COO-和CO-O之间共享共振结构的等同氧原子. 这些原子是羰基氧, 与蛋白质骨架中的氧相似, 但实际上与去质子化天冬氨酸侧链中的氧能更好地匹配. 在这种情况下, 我们将这些原子归类为O2类型.

3.2.2) 定义一个不相似的氧原子的原子类型

氧杂蒽O15中的氧原子是芳香性的, 因此与FF99SB中已经存在的任何氧原子类型都不相符, 因为它们都是SP2型羰基或SP3型醚氧. 为一个新的原子类型生成参数超出了本教程的范围, 幸好此特定原子类型的参数已被发布(VanBeek 等人, 生物物理学. J. 2007, 92, 4168-4178 PDF supp mat). 所以我们可以直接利用这些参数. 在这时, 我们需要考虑的是在当前的力场文件中选择正确的原子类型. 按照VanBeek文献中的术语, 我们将它定义为OA.

剩下的原子类型选择如下图所示.

图3.5

3.3) 保存新的库文件

现在我们已经输入了xleap所需的所有参数来充分识别FAM5残基, 我们可以保存一个库文件, 使得xleap在以后也能识别这个残基. 这非常重要因为我们不必每次重复上述所有步骤将FAM5导入到xleap程序. 为此, 输入以下命令以保存一个新的pdb.

saveoff FAM5 fam5.lib
savepdb FAM5 fam5_leap.pdb

3.4) 创建一个frcmod文件.

需要一个frcmod文件来为我们添加的新的OA原子类型定义质量和范德华参数, 并提供所有在标准FF99SB力场中不存在的键, 角和二面角参数. 我们将使用leap程序来帮助我们识别缺失的参数. 这将需要几个步骤, 因为leap检查功能只能识别缺失的键和角参数, 并且寻找丢失的二面角, 在我们提供缺失的键和角参数之后, 我们将需要尝试保存prmtop和inpcrd文件. 如果你已经跳过了, 你应关闭它(确保你保存了你的lib和pdb文件). 由于我们不需要图形界面, 所以我们可以使用更快, 更方便的tleap.

$AMBERHOME/exe/tleap -s -f $AMBERHOME/dat/leap/cmd/ oldff/leaprc.ff99SB

然后我们可以加载我们上面创建的lib文件:

>loadoff fam5.lib

现在我们可以检查fam5单元:

>check FAM5

Leap将显示缺失的键和角度参数:

Checking 'FAM5'....
ERROR: The unperturbed charge of the unit: -2.307860 is not integral.
WARNING: The unperturbed charge of the unit: -2.307860 is not zero.
ERROR: The perturbed charge: -2.307860 is not integral.
WARNING: The perturbed charge: -2.307860 is not zero.
Checking parameters for unit 'FAM5'.
Checking for bond parameters.
Could not find bond parameter for: CA - O
No bond parameter for: CA - O
Could not find bond parameter for: CA - OA
No bond parameter for: CA - OA
Could not find bond parameter for: OA - CA
No bond parameter for: OA - CA
Could not find bond parameter for: CA - O
No bond parameter for: CA - O
Checking for angle parameters.
Could not find angle parameter: O - C - CA
Can't find angle parameter: O - C - CA
Could not find angle parameter: CA - C - N
Can't find angle parameter: CA - C - N
Could not find angle parameter: CA - C - O2
Can't find angle parameter: CA - C - O2
Could not find angle parameter: CA - C - O2
Can't find angle parameter: CA - C - O2
Could not find angle parameter: CA - CA - OA
Can't find angle parameter: CA - CA - OA
Could not find angle parameter: CA - CA - O
Can't find angle parameter: CA - CA - O
Could not find angle parameter: O - CA - CA
Can't find angle parameter: O - CA - CA
Could not find angle parameter: CA - CA - OA
Can't find angle parameter: CA - CA - OA
Could not find angle parameter: CA - OA - CA
Can't find angle parameter: CA - OA - CA
Could not find angle parameter: OA - CA - CA
Can't find angle parameter: OA - CA - CA
Could not find angle parameter: OA - CA - CA
Can't find angle parameter: OA - CA - CA
Could not find angle parameter: CA - CA - O
Can't find angle parameter: CA - CA - O
Could not find angle parameter: O - CA - CA
Can't find angle parameter: O - CA - CA
There are missing parameters.
check: Errors: 2 Warnings: 2

前两个错误和警告是因为我们的FAM5单位没有整数电荷. 正如前面所讨论的那样, 这是我们所期望的, 因为AMBER中的5’残基不应该具有整数电荷. 下面的列表是所有我们缺少的参数, 虽然有相当多的冗余. 例如, 提示中的最后两个参数是CA-CA-O和O-CA-CA指的是相同的角度, 只是在原子的顺序不同. 所以, 去掉这些冗余剩下的参数列表就是我们需要定义的:

CA - O
CA - OA
CA - C - O(从O - C - CA重新排序)
CA - C - N
CA - C - O2
CA - CA - OA
CA - CA - O
CA - OA - CA

为这些参数提供值是非常重要的. 如果幸运的话, 在工作环境中这些缺失的参数有与现有参数非常相似的原子, 我们可以通过”类似”来定义它们. 例如, 在广义AMBER力场中存在一个ca-ca-o参数, 我们可以将它用到CA-CA-O参数. 这种参数分配必须小心, 彻底的分配超出了本教程的范围. 从前面提到的VanBeek文献找相关参数(VanBeek 等生物物理学J. 2007, 92, 4168-4178 PDF SUPP mat), 我们提出了在frcmod文件中的如下第一关, fam5_incomplete.frcmod.

VanBeek等人 Biophys J.(2007)92,4168-4178
质量
OA  16.00  0.465, 与gaff 的os和parm99相同
键
CA-OA  372.40  1.373与gaff 的ca- os相同
O  648.00  1.214与gaff 的c -o相同
角
CA-CA-OA 69.800  119.200与gaff的 ca-ca-os相同
CA-CA-O  70.900  123.430与gaff的ca-ca-o 相同
CA-C -0  68.700  123.440与gaff的 ca-c -o 相同
CA-C -O2  68.700  123.440与gaff的g-ca-c -o相同
CA-OA-CA  63.600  118.960与GAFF的 ca-os-ca 相同
CA-C -N 69.400 112.03与GAFF的 ca-c -n相同
二面角
异二面角
非键信息
OA  1.6837  0.1700  OPLS醚等与gaff的 os相同

注意这些参数提供了分子的合理描述, 但是如果需要高度准确表示荧光素的话, 则需要对这些参数进行优化. 还要注意, 我们还没有包含任何二面角参数, 因为我们还没有算出这些缺失的参数. 最后, 即使check命令没有提供有关缺少原子类型OA的警告消息, 我们也知道这个问题, 所以我们也引入了这个原子类型的质量和非键参数.

我们现在加载fam5_incomplete.frcmod文件, 然后找出需要定义的二面角参数:

> loadamberparams fam5_incomplete.frcmod

不幸的是, 检查命令不检查二面角, 所以我们只是先试着保存prmtop和inpcrd文件, 然后leap就会列出缺少的二面角参数.

> saveamberparm FAM5 prmtop inpcrd

Leap将列出缺失的键和角度参数:

 ** No torsion terms for  CA-CA-OA-CA
 ** No torsion terms for  CA-CA-OA-CA
 ** No torsion terms for  CA-OA-CA-CA
 ** No torsion terms for  CA-OA-CA-CA

在这种情况实际上只缺失一个二面角参数(CA-CA-OA-CA), 因为上述四个都描述的是相同的二面角. 所以我们只需要在frcmod文件的DIHE部分添加一个条目:

CA-CA-OA-CA  4  14.500  180.000  2.000与gaff的 X-ca-ca-X相同, 因为M.Zwier在手动修改时发现~18

使用文本编辑器更新fam5_incomplete.frcmod文件并将其保存为fam5.frcmod. 现在将其加载到Leap中, 并再次尝试saveamberparm命令作为最终检查; 实际上我们没有使用生成的prmtop和inpcrd文件.

我们已经成功保存了我们的FAM5库文件并创建了相关的frcmod文件, 现在可以将其附加到我们的polyAT十聚体中.

第4部分 将FAM5分子连接到polyAT十聚体的5’末端

在第4节中, 我们将FAM5分子附加到我们的polyAT十聚体的5’末端(polyAT.pdb). 我们将使用与第2节中描述的类似的过程, 其中涉及广泛使用Sirius可视化程序. 下面描述了FAM分子与linker的连接(图4.1). 然而, 将分子连接到任何DNA或RNA双链体的5’末端存在几个问题. 第一个问题是polyAT十聚体的5’末端丢失的磷酸基团. 为了连接我们的FAM分子, 我们必须通过将末端腺嘌呤变成内部腺嘌呤来将缺失的磷酸酯添加到十聚体中. 我们将使用xleap来修改polyAT的5’末端, 因为它可以自动修改腺嘌呤残基. 第二个问题是FAM和polyAT分子之间缺乏重复的原子(因为我们必须删除linker上的磷酸基团), 这可能会被用来叠加键. 因此, 我们必须手动将linker移动到相对于polyAT十聚体来说相对正确键的位置. 须注意的是, 任何手动构建的键都可能会变形, 因此在跑生产相的MD之前需要小心地将结构能量最小化.

图4.1

4.1) 准备polyAT十聚体

如前所述, 我们的第一步是将磷酸基团添加到polyAT十聚体的5’末端. Xleap通过使用它的”loadpdbuseseq”命令来简化这个过程. 这个命令使我们能够定义我们想要加载到程序中的确切序列, 并指定一个pdb文件来对应这个序列. 让我们开始执行xleap程序.

$AMBERHOME/exe/xleap &

在xleap窗口中, 我们将加载ff99SB力场, 以便识别我们的DNA.

source oldff/leaprc.ff99SB

接下来, 我们将定义我们想要加载到xleap的序列. 这是最重要的一步, 因为我们将使用此过程来指定我们的残基修饰. 而不是指定在pdb中的DA5残基, 我们将告诉xleap加载这个残基的DA模板. 这是一个内部残基, 所以xleap将添加我们连接linker所需的额外磷酸脂.

polyATseq = {DA DA DA DA DA DA DA DA3 DT5 DT DT DT DT DT DT DT DT3}

通过指定上述顺序时, 我们将告诉xleap加载我们的polyAT.pdb文件, 但忽略pdb中的残基名称, 而是使用我们提供的序列.

polyAT = loadpdbuseseq polyAT.pdb polyATseq

以下是xleap窗口的输出内容:

载入PDB文件: ./polyAT.pdb使用序列polyATseq 匹配pdb残基 - >序列模板 res pdb template 1* A5 DA 2 A DA 3 A DA 4 A DA 5 A DA 6 A DA 7 A DA 8 A DA 9 A DA 10 A3 DA3 11 T5 DT5 12 T DT 13 T DT 14 T DT 15 T DT 16 T DT 17 T DT 18 T DT 19 T DT 20 T3 DT3

• =可能不匹配; 共1 (即pdb名称不是模板的子字符串) 在残基内创建了一个名为H5T的新原子: .R <DA 1> 添加缺失的重原子: .R <DA 1> .A <P 1> 添加缺失的重原子: .R <DA 1> .A <O1P 2> 在文件中添加缺失的重原子: .R <DA 1> .A <O2P 3> 总原子: 438 根据残基模板, 增加了203个失去的原子: 3个重 200H /孤对电子 文件包含1个不在残基模板中的原子

发生了什么？

Xleap将对比我们指定的序列与pdb文件的序列, 并指出pdb文件中的5’腺嘌呤(A5)与序列中指定的内部腺嘌呤(DA)之间可能不匹配的地方. 这正是我们所期望的. Xleap还使用重原子P和O将磷酸酯添加到我们的分子中. 然而, 我们应该注意这个新的原子”H5T”. 这意味着xleap不能识别. 这个原子是在原来的pdb文件中一个额外的质子, 且与内部腺嘌呤中的原子不对应. 为了检查这个原子, 我们使用编辑窗口来查看这个结构.

edit polyAT

乍一看, 我们注意到在polyAT十聚体的5’末端有一个非键合原子(在下图中圈出). 要显示原子名称, 请选择显示>名称. 你也可以放大分子以获得更好的观察. 该分子的原子名称H5T与由xleap创建的原子相匹配. Xleap聪明地将缺少的原子添加到已指定的残基, 但是, 它不会删除pdb文件中存在的额外原子. H5T原子是5’腺嘌呤残基多余的, 应该删除. 你可以使用位于工具栏上的橡皮擦来删除此原子.

图4.2

现在我们已经将磷酸基团添加到了我们的polyAT十聚体的5’腺苷酸末端, 我们可以保存一个新的pdb文件并附加我们的FAM分子. 关闭编辑窗口并保存polyAT的新的pdb文件polyAT_edit.pdb.

savepdb polyAT polyAT_edit.pdb

注意: 这里的另一种方法是在将其加载到leap之前把H5T原子从polyAT.pdb文件中删除掉.

4.2) 使用Sirius将FAM和polyAT分子合并在一起

在合并polyAT和FAM5的pdb文件时, 我们需要确保这两个单元相对于对方位置合理. 就像在第二部分一样, 我们将再次使用Sirius. 启动Sirius应用程序并打开polyAT_edit.pdb和fam5_leap.pdb. (注意, 如果在图形窗口中没有看到结构, 则可能需要单击雾化按钮并打开和关闭深度提示选择使其显示). 这些结构将会彼此重叠. 我们需要通过移动其中一个结构, 同时保持另一个不动. 原则上, 只要我们在整个过程中都保持一致, 那么我们移动哪个都不重要; 我们将选择移动fam5部分. 只要在Sirius中移动fam5部分, 单击Set Separate Motion的按钮.

在打开的对话框中, 选择fam5_leap, 然后单击确定. 现在在Sirius平移和旋转只影响fam5分子. 你应注意FAM5的位置, 使氧原子与磷酸酯如下图4.3和4.4所示紧紧靠近. 在旋转FAM 5单位(鼠标左键), 平移FAM 5(XY右键和Z中心键)之间往往需要重复, 然后旋转整个体系以获得更好的视图 - 通过在选择”set separate motion button “之后立即选择”set combined motion “. 要注意保持PO键合理的键长, 避免空间碰撞. 这需要一定的耐心. 注意使用正确的DNA残基去结合.

图4.3

图4.4

一旦你对这个视图的角度感到满意, 请选择”文件” - >”保存结构”. Sirius会提示你想保存哪个对象. 选择FAM5_leap并点击保存. 当提示保存新的坐标说是的. 然后把它保存为fam5_newcoord.pdb.

不好的是, 由Sirius生成的pdb文件具有与Leap不兼容的连接信息, 所以我们需要删除它. 使用文本编辑器去除连接信息和”END”标签, 并将其保存为fam5_m.pdb. 接下来, 我们需要将fam5和polyAT pdb文件合并在一起polyAT_edit.pdb fam5_polyAT_edit.pdb:

cat fam5_m.pdb polyAT_edit.pdb > fam5_polyAT_edit.pdb

4.3) 将组合的pdb加载到Leap并创建prmtop和inpcrd文件

为了读取组合的pdb, 我们需要做一些事情. 首先, 我们需要加载DNA的参数以及我们在第3节中创建的fam5的frcmod文件. 我们还需要加载也在第3节中创建的fam5的lib文件. 但是, 我们需要稍作处理, 我们的fam5库文件是以FAM5作为残留名称, 而pdb文件是FAM(这是因为Sirius只支持3个字符残基名称). 我们将使用leap的复制命令来解决这个问题(或者你可以编辑pdb并用’FAM5’替换’FAM’).

xleap -f -s $AMBERHOME/dat/leap/cmd/oldff/leaprc.ff99SB
>loadoff fam5.lib
>loadamberparams fam5.frcmod
>FAM = copy FAM5
>mol = loadpdb fam5_polyAT_edit.pdb

这是leap将会提示一些警告, 如下所示:

加载PDB文件: ./fam5_polyAT_edit.pdb
警告: 在pdb文件中的1号残基更改名称;
这个残基被分成FAM和DA.
有1个残基有命名警告.
有分裂的残基;
残留序列号将不对应于pdb中的序列号.
单方面的连接. 残基: 缺少连接1原子.
文件中的总原子数: 697

我们实际上可以忽略这些警告. 第一个是因为FAM和第一个DA残差都有1的残基ID, 因为我们只是简单的附加文件. 然而, Leap已经检测到这个冲突并更新了编号. 第二个警告是FAM残基不含有尾部原子, 因此没有原子与第一个DA残基的头部原子结合. 我们可以在第3节中通过定义一个尾原子来解决这个问题. 但是, 我们可以简单通过一个手动键合的命令来解决这个问题:

>bond mol.1.57 mol.2.1

现在你可以用sander或pmemd来保存气相/隐式溶剂的prmtop和inpcrd文件:

>saveamberparm mol fam5_polyAT_vac.prmtop fam5_polyAT_vac.inpcrd

最后如果你想跑显示溶剂模拟, 你可以添加反离子和溶剂, 就像在前面的教程中所做的那样.