今天开始学习跟着官网的tutorial学习Gromacs2020.3的操作和运行。tutorial本尊|
Gromacs的安装很简单，我是在conda下面安装的的，直接就是2020.3的版本。
Tutorial 1：Molecular Modeling Practical 分子建模
模型分子 泛素连接酶E2
1.简单背景 👇
泛素连接酶也叫E2酶，参与泛素途径的第二步骤。这一反应属于蛋白的转录后修饰，实在蛋白上面共价结合1个或者数个泛素单体。被泛素标记的蛋白会被蛋白酶降解，重利用。另外泛素化本身也在很多蛋白的结构稳定，功能的调整以及包内转运方面发挥重要作用。
这个泛素化的过程是需要ATP供能的，至少有三个酶的参与，泛素活化酶E1，泛素结合酶E2以及泛素连接酶E3。这里面参与的蛋白发生突变会影响整个通路，产生诸如肿瘤、帕金森氏和一些西南哦血管的疾病。泛素化标签的传递主要特异性的以来E2-E3直接的相互作用。
人类的基因组存在37中功能不同的E2结合酶和700余种E3连接酶。E3的多样性就要求每个E2会链接不同构形的E3连接酶。泛素结合酶也就是E2大家族，存在一个高度保守的泛素结合域（UBC）。E3就是通过这个保守区域和E2结合的，UBC包括helix1，loop1和loop2。本教程在这里研究野生型E2的分子动力学。
2.Molecular Dynamics ⚙
经典的分子动力学模拟使用牛顿运动方程计算粒子的轨迹。计算系统里的每一个粒子和其他粒子的相互作用，计算每个粒子在这个力场中所获得的加速度及方向，用这些数据来决定下一步的位置。这样基于实验数据的分子动力学模拟的结果对于理解分子对接原理和新的假说非常有帮助。但是，受限于计算能力，能模拟的体系大小和时间都是很有限的。
3.Gromacs ☯
本教程将用到Gromacs2020来找事E2的分子动力学模拟的过程。Gromacs是一套基于经典力学进行分子模拟的免费软件，运行在linux command-line下
4.开始
4.1 初始结构
我们需要下载结构，PDBID:1Y6L 和 3BZH 可以使用Pymol或者其他任何可视化软件进行观察。
4.2 准备工作
A simplified version of the 1Y6L pdb file is waiting for you as E2Bwt.pdb.在做之前有必要说一下命名规则，每次做完的操作得到的结果都要体现在结果文件名上，这样有利于结果的整理个查找。为了这个教程对你的蛋白结果也有通用性，请先把你的蛋白也改名为 protein.pdb 并转移到新的文件夹下，确保这个文件夹下只有这一个文件。
4.3 文件结构转换和拓扑文件的生成
gmx pdb2gmx -f protein.pdb -o protein.gro -p protein.top -ignh
选择合适的力场文件:这个例子来我选的是 53a6
Select the Force Field:
From '/home/lxb/miniconda3/envs/MD/share/gromacs/top':
1: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
2: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
3: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)
4: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
5: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)
6: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)
7: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
8: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)
9: GROMOS96 43a1 force field
10: GROMOS96 43a2 force field (improved alkane dihedrals)
11: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
12: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
13: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
14: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40,, 843-856, DOI: 10.1007/s00249-011-0700-9)
15: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
选择溶液模型模式,这里面我们选的是SPC简单点电荷的模式。也可以选择场模式，这个在以后的贴里面再讨论。
Select the Water Model:
1: SPC simple point charge, recommended
2: SPC/E extended simple point charge
3: None
瞬间结束，我们得到三个文件
├── posre.itp # 另一个结果文件
├── protein.gro # gro文件
├── protein.pdb
└── protein.top # 拓扑结构文件
You can always generate a pdb file from a gro file using the following command: editconf -f XXXXX.gro -o XXXXX.pdb
4.3 蛋白结构的能量最小化
现在结构已经用正确的格式展示出来了。但是，这个结构里面可能存在某些不合理的地方，例如确实或者多氢的情况，这个是很难避免的，或者晶体衍射数据整合的时候原子在某些位置过于集中这样的不问他区域，因此我们要做一个能量最小化的操作，把这些局部的过高能量释放掉。我们用minim.mdp
; Lines starting with ';' ARE COMMENTS
; Everything following ';' is also comment
title = Energy Minimization ; Title of run
; The following line tell the program the standard locations where to find certain files
cpp = /lib/cpp ; Preprocessor
; Define can be used to control processes
define = -DFLEXIBLE
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 10.0 ; Stop minimization when the maximum force < 1.0 kJ/mol
nsteps = 5000 ; Maximum number of (minimization) steps to perform
nstenergy = 1 ; Write energies to disk every nstenergy steps
energygrps = System ; Which energy group(s) to write to disk
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = Reaction-Field ; Treatment of long range electrostatic interactions
epsilon_rf = 78
rcoulomb = 1.4 ; long range electrostatic cut-off
rvdw = 1.4 ; long range Van der Waals cut-off
constraints = none ; Bond types to replace by constraints
pbc = xyz ; Periodic Boundary Conditions (yes/no)
来进行这一5000steps的能量最小化操作。
grompp -v -f minim.mdp -c protein.gro -p protein.top -o protein-EM-vacuum.tpr
这一步之后会生成tpr文件，这一文件包括所有参数，蛋白格点文件信息。我们直接用mdrun这个文件投任务。
gmx mdrun -v -deffnm protein-EM-vacuum -c protein-EM-vacuum.gro
因为是要做能量最小化，因此不需要过程文件，我们只需要能量最小的终点结构即可。
4.4 周期性边界
目前这个系统已经做到了能量最小化，但是目前这个系统还是真空的还没有空间和任何溶剂的信息。所以我们需要把这个蛋白放到一个特定的空间里面，并加上溶剂。为了在有限空间模拟更多的功能，我们把边界周期化。用editconf进行。通常情况下我们采用菱形十二面体因为他近似球体。所以蛋白可以在这里面自由旋转。
gmx editconf -f protein-EM-vacuum.gro -o protein-PBC.gro -bt dodecahedron -d 1.2
4.5 加入水溶液
现在我们可以在这个空间里面充满水溶液，有两种方法，一个是按照水的密度直接充满整个空间。另一个是用水把蛋白包裹起来。这样可以节省计算资源，但是这样的话可能会要求空间计算的要少紧密一些否则很难保证谁的密度。我们采用和水密度相同的方式充满整个空间
gmx solvate -cp protein-PBC.gro -cs spc216.gro -p protein.top -o protein-water.gro
结果中
Output configuration contains 32889 atoms in 10613 residues
Volume : 340.649 (nm^3)
Density : 1001.2 (g/l)
Number of solvent molecules: 10465
水的密度是1吧 齐活。
加好了水 但是 不一定电荷是平的，有可能有一两个正电子 或者 负电，我们用NaCl把电荷补平 这样才科学。
grompp -v -f minim.mdp -c protein-water.gro -p protein.top -o protein-water.tpr
看到没？ 有两个正电荷需要补齐。
NOTE 2 [file protein.top, line 9505]:
System has non-zero total charge: 2.000000
Total charge should normally be an integer. See
http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
for discussion on how close it should be to an integer.
要用NaCl来代替水分子。毕竟是溶于水的嘛~
genion -s protein-water.tpr -o protein-solvated.gro -conc 0.15 -neutral -pname NA+ -nname CL-
……
eading file protein-water.tpr, VERSION 2020.3 (single precision)
Reading file protein-water.tpr, VERSION 2020.3 (single precision)
Will try to add 31 NA+ ions and 33 CL- ions.
Select a continuous group of solvent molecules
Group 0 ( System) has 32889 elements
Group 1 ( Protein) has 1494 elements
Group 2 ( Protein-H) has 1163 elements
Group 3 ( C-alpha) has 148 elements
Group 4 ( Backbone) has 444 elements
Group 5 ( MainChain) has 593 elements
Group 6 ( MainChain+Cb) has 733 elements
Group 7 ( MainChain+H) has 730 elements
Group 8 ( SideChain) has 764 elements
Group 9 ( SideChain-H) has 570 elements
Group 10 ( Prot-Masses) has 1494 elements
Group 11 ( non-Protein) has 31395 elements
Group 12 ( Water) has 31395 elements
修改top文件的最后几行，把NA CL带上把相应水分子去掉。
[ molecules ]
; Compound #mols
Protein_chain_A 1
SOL 10401
NA 31
CL 33
带着水分子继续 做能量最小化，舒展下水分子筋骨
gmx grompp -v -f minim.mdp -c protein-solvated.gro -p protein.top -o protein-EM-solvated.tpr
mdrun -v -deffnm protein-EM-solvated -c protein-EM-solvated.gro
4.5 生成tpr跑MD
我们需要一个参数文件 pr.mdp
gmx grompp -v -f pr.mdp -c protein-EM-solvated.gro -r protein-EM-solvated.gro -p protein.top -o protein-PR.tpr
来进行分子动力学的模拟计算。
投递这个任务可以用nohup 或者集群的任意投递方式即可
mdrun -v -deffnm protein-PR
网友评论