仓库源文，站点原文

layout: post title: GROMACS非标准残基教程2：芋螺毒素小肽实例 categories:

科 tags:
gmx
amber chem: true

2017-09-22 22:23:42 整理: 李睿; 补充: 康文渊; 修订: 李继存

在本教程的第一部分中, 我们整理了GROMACS处理非标准残基的理论流程, 但如果你不熟悉GROMACS的力场结构, 单单根据那个流程来处理非标准残基仍然存在很大困难. 为此, 在本教程的第二部分中, 我们以一个具体实例来说明使用GROMACS时, 如何处理非标准残基, 如何在GROMACS力场文件中添加新的残基数据.

我们要处理的蛋白是GIC, 它是一种具有生理活性的芋螺毒素小肽, 其主要作用是抑制乙酰胆碱受体的开放, 因此有很大的潜力被开发成为一种止痛药. 在开发过程中, 有文章报道使用伽马氨基丁酸修饰后GIC会改变其活性. 这种修饰后的小肽被称为GID, 其中的伽马氨基丁酸残基CGU是一种非标准氨基酸残基. 因此, 如果想要使用GROMACS对GID进行分子动力学模拟, 就必须处理非标准残基CGU.

由于要模拟的体系是蛋白, 因此我们选用amber99sb-ildn力场. 但教程中的处理方法同样适用于其他力场.

第一步获取蛋白的PDB文件

GID的PDB登录号为1MTQ, 下载其PDB文件1mtq.pdb.

下载后, 按常规步骤查看并清洗pdb文件, 去除杂原子, 去除重复模型, 只保留一个蛋白模型, 将其保存为gid.pdb. GID结构如下

<figure><script>var Mol1=new ChemDoodle.TransformCanvas3D('Mol-1',642,396);Mol1.specs.shapes_color = '#fff';Mol1.specs.backgroundColor = 'black';Mol1.specs.compass_display = true;Mol1.specs.set3DRepresentation('Ball and Stick');Mol1.specs.projectionPerspective_3D = false;Mol1.specs.proteins_ribbonCartoonize = true;Mol1.handle = null;Mol1.timeout = 15;Mol1.specs.crystals_unitCellLineWidth = 1.5;Mol1.specs.nucleics_residueColor = 'rainbow';Mol1.specs.proteins_residueColor= 'rainbow';Mol1.startAnimation = ChemDoodle._AnimatorCanvas.prototype.startAnimation;Mol1.stopAnimation = ChemDoodle._AnimatorCanvas.prototype.stopAnimation;Mol1.isRunning = ChemDoodle._AnimatorCanvas.prototype.isRunning;Mol1.dblclick = ChemDoodle.RotatorCanvas.prototype.dblclick;Mol1.nextFrame = function(delta){var matrix = [];ChemDoodle.lib.mat4.identity(matrix);var change = delta*Math.PI/15000;ChemDoodle.lib.mat4.rotate(matrix,change,[1,0,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,1,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,0,1]);ChemDoodle.lib.mat4.multiply(this.rotationMatrix, matrix)};var Fmol='HEADER TOXIN 22-SEP-02 1MTQ \nTITLE THREE-DIMENSIONAL SOLUTION STRUCTURE OF ALPHA-CONOTOXIN GID \nTITLE 2 BY NMR SPECTROSCOPY \nCOMPND MOL_ID: 1; \nCOMPND 2 MOLECULE: ALPHA-CONOTOXIN GID; \nCOMPND 3 CHAIN: A; \nCOMPND 4 ENGINEERED: YES \nSOURCE MOL_ID: 1; \nSOURCE 2 SYNTHETIC: YES; \nSOURCE 3 OTHER_DETAILS: THIS SEQUENCE OCCURS IN NATURALLY IN CONUS \nSOURCE 4 GEOGRAPHUS VENOM. \nKEYWDS ALPHA-HELIX, TOXIN \nEXPDTA SOLUTION NMR \nNUMMDL 20 \nAUTHOR A.NICKE,M.L.LOUGHNAN,E.L.MILLARD,P.F.ALEWOOD,D.J.ADAMS, \nAUTHOR 2 N.L.DALY,D.J.CRAIK,R.J.LEWIS \nREVDAT 2 24-FEB-09 1MTQ 1 VERSN \nREVDAT 1 11-FEB-03 1MTQ 0 \nJRNL AUTH A.NICKE,M.L.LOUGHNAN,E.L.MILLARD,P.F.ALEWOOD, \nJRNL AUTH 2 D.J.ADAMS,N.L.DALY,D.J.CRAIK,R.J.LEWIS \nJRNL TITL ISOLATION, STRUCTURE, AND ACTIVITY OF GID, A NOVEL \nJRNL TITL 2 ALPHA 4/7-CONOTOXIN WITH AN EXTENDED N-TERMINAL \nJRNL TITL 3 SEQUENCE \nJRNL REF J.BIOL.CHEM. V. 278 3137 2003 \nJRNL REFN ISSN 0021-9258 \nJRNL PMID 12419800 \nJRNL DOI 10.1074/JBC.M210280200 \nREMARK 1 \nREMARK 2 \nREMARK 2 RESOLUTION. NOT APPLICABLE. \nREMARK 3 \nREMARK 3 REFINEMENT. \nREMARK 3 PROGRAM : CNS 1.0 \nREMARK 3 AUTHORS : BRUNGER ET AL. \nREMARK 3 \nREMARK 3 OTHER REFINEMENT REMARKS: THE STRUCTURES ARE BASED ON A TOTAL \nREMARK 3 OF 183 NOE DERIVED DISTANCE RESTRAINTS, 15 DIHEDRAL ANGLE \nREMARK 3 RESTRAINTS. \nREMARK 4 \nREMARK 4 1MTQ COMPLIES WITH FORMAT V. 3.15, 01-DEC-08 \nREMARK 100 \nREMARK 100 THIS ENTRY HAS BEEN PROCESSED BY RCSB ON 27-SEP-02. \nREMARK 100 THE RCSB ID CODE IS RCSB017175. \nREMARK 210 \nREMARK 210 EXPERIMENTAL DETAILS \nREMARK 210 EXPERIMENT TYPE : NMR \nREMARK 210 TEMPERATURE (KELVIN) : 287 \nREMARK 210 PH : 2.8 \nREMARK 210 IONIC STRENGTH : 0 \nREMARK 210 PRESSURE : AMBIENT \nREMARK 210 SAMPLE CONTENTS : 2.5MM PEPTIDE; 90% H2O, 10% \nREMARK 210 D2O \nREMARK 210 \nREMARK 210 NMR EXPERIMENTS CONDUCTED : 2D TOCSY, 2D NOESY, DQF-COSY, \nREMARK 210 E-COSY \nREMARK 210 SPECTROMETER FIELD STRENGTH : 750 MHZ \nREMARK 210 SPECTROMETER MODEL : DRX \nREMARK 210 SPECTROMETER MANUFACTURER : BRUKER \nREMARK 210 \nREMARK 210 STRUCTURE DETERMINATION. \nREMARK 210 SOFTWARE USED : XWINNMR 2.6, XEASY 1.3.7, \nREMARK 210 DYANA 1.5 \nREMARK 210 METHOD USED : INITIAL STRUCTURES WERE \nREMARK 210 CALCULATED USING TORSION ANGLE \nREMARK 210 DYNAMICS AND STRUCTURES WERE \nREMARK 210 REFINED IN A WATER SHELL USING \nREMARK 210 CARTESIAN DYNAMICS. \nREMARK 210 \nREMARK 210 CONFORMERS, NUMBER CALCULATED : 50 \nREMARK 210 CONFORMERS, NUMBER SUBMITTED : 20 \nREMARK 210 CONFORMERS, SELECTION CRITERIA : ALL CALCULATED STRUCTURES \nREMARK 210 SUBMITTED, STRUCTURES WITH THE \nREMARK 210 LOWEST ENERGY \nREMARK 210 \nREMARK 210 BEST REPRESENTATIVE CONFORMER IN THIS ENSEMBLE : 16 \nREMARK 210 \nREMARK 210 REMARK: THIS STRUCTURE WAS DETERMINED USING STANDARD 2D \nREMARK 210 HOMONUCLEAR TECHNIQUES. \nREMARK 215 \nREMARK 215 NMR STUDY \nREMARK 215 THE COORDINATES IN THIS ENTRY WERE GENERATED FROM SOLUTION \nREMARK 215 NMR DATA. PROTEIN DATA BANK CONVENTIONS REQUIRE THAT \nREMARK 215 CRYST1 AND SCALE RECORDS BE INCLUDED, BUT THE VALUES ON \nREMARK 215 THESE RECORDS ARE MEANINGLESS. \nREMARK 500 \nREMARK 500 GEOMETRY AND STEREOCHEMISTRY \nREMARK 500 SUBTOPIC: TORSION ANGLES \nREMARK 500 \nREMARK 500 TORSION ANGLES OUTSIDE THE EXPECTED RAMACHANDRAN REGIONS: \nREMARK 500 (M=MODEL NUMBER; RES=RESIDUE NAME; C=CHAIN IDENTIFIER; \nREMARK 500 SSEQ=SEQUENCE NUMBER; I=INSERTION CODE). \nREMARK 500 \nREMARK 500 STANDARD TABLE: \nREMARK 500 FORMAT:(10X,I3,1X,A3,1X,A1,I4,A1,4X,F7.2,3X,F7.2) \nREMARK 500 \nREMARK 500 EXPECTED VALUES: GJ KLEYWEGT AND TA JONES (1996). PHI/PSI- \nREMARK 500 CHOLOGY: RAMACHANDRAN REVISITED. STRUCTURE 4, 1395 - 1400 \nREMARK 500 \nREMARK 500 M RES CSSEQI PSI PHI \nREMARK 500 3 ARG A 2 34.82 -82.74 \nREMARK 500 5 SER A 7 36.97 -83.90 \nREMARK 500 6 ARG A 2 -42.10 -170.65 \nREMARK 500 9 ARG A 2 -76.92 -90.42 \nREMARK 500 9 SER A 7 39.08 -88.91 \nREMARK 500 10 SER A 7 30.23 -86.57 \nREMARK 500 11 SER A 7 44.93 -84.67 \nREMARK 500 14 ARG A 2 -36.14 77.91 \nREMARK 500 14 SER A 7 41.23 -80.22 \nREMARK 500 15 ARG A 2 -10.69 -161.81 \nREMARK 500 15 SER A 7 30.57 -95.57 \nREMARK 500 20 SER A 7 44.90 -87.46 \nREMARK 500 \nREMARK 500 REMARK: NULL \nREMARK 999 \nREMARK 999 AN APPROPRIATE SEQUENCE DATABASE MATCH WAS NOT \nREMARK 999 FOUND USING BLAST. \nSEQRES 1 A 19 ILE ARG ASP CGU CYS CYS SER ASN PRO ALA CYS ARG VAL \nSEQRES 2 A 19 ASN ASN HYP HIS VAL CYS \nMODRES 1MTQ CGU A 4 GLU GAMMA-CARBOXY-GLUTAMIC ACID \nMODRES 1MTQ HYP A 16 PRO 4-HYDROXYPROLINE \nHET CGU A 4 17 \nHET HYP A 16 15 \nHETNAM CGU GAMMA-CARBOXY-GLUTAMIC ACID \nHETNAM HYP 4-HYDROXYPROLINE \nHETSYN HYP HYDROXYPROLINE \nFORMUL 1 CGU C6 H9 N O6 \nFORMUL 1 HYP C5 H9 N O3 \nHELIX 1 1 ASN A 8 ASN A 15 1 8 \nSSBOND 1 CYS A 5 CYS A 11 1555 1555 2.03 \nSSBOND 2 CYS A 6 CYS A 19 1555 1555 2.03 \nLINK C ASP A 3 N CGU A 4 1555 1555 1.32 \nLINK C CGU A 4 N CYS A 5 1555 1555 1.33 \nLINK C ASN A 15 N HYP A 16 1555 1555 1.34 \nLINK C HYP A 16 N HIS A 17 1555 1555 1.33 \nCRYST1 1.000 1.000 1.000 90.00 90.00 90.00 P 1 1 \nORIGX1 1.000000 0.000000 0.000000 0.00000 \nORIGX2 0.000000 1.000000 0.000000 0.00000 \nORIGX3 0.000000 0.000000 1.000000 0.00000 \nSCALE1 1.000000 0.000000 0.000000 0.00000 \nSCALE2 0.000000 1.000000 0.000000 0.00000 \nSCALE3 0.000000 0.000000 1.000000 0.00000 \nMODEL 1 \nATOM 1 N ILE A 1 -14.555 7.084 -3.511 1.00 0.00 N \nATOM 2 CA ILE A 1 -13.512 6.116 -3.816 1.00 0.00 C \nATOM 3 C ILE A 1 -12.411 6.178 -2.763 1.00 0.00 C \nATOM 4 O ILE A 1 -11.919 7.258 -2.432 1.00 0.00 O \nATOM 5 CB ILE A 1 -12.902 6.372 -5.214 1.00 0.00 C \nATOM 6 CG1 ILE A 1 -14.003 6.381 -6.279 1.00 0.00 C \nATOM 7 CG2 ILE A 1 -11.851 5.317 -5.542 1.00 0.00 C \nATOM 8 CD1 ILE A 1 -13.508 6.731 -7.667 1.00 0.00 C \nATOM 9 H ILE A 1 -14.331 7.862 -2.955 1.00 0.00 H \nATOM 10 HA ILE A 1 -13.954 5.131 -3.809 1.00 0.00 H \nATOM 11 HB ILE A 1 -12.416 7.335 -5.200 1.00 0.00 H \nATOM 12 HG12 ILE A 1 -14.456 5.403 -6.327 1.00 0.00 H \nATOM 13 HG13 ILE A 1 -14.754 7.108 -6.002 1.00 0.00 H \nATOM 14 HG21 ILE A 1 -11.941 5.030 -6.580 1.00 0.00 H \nATOM 15 HG22 ILE A 1 -12.003 4.450 -4.915 1.00 0.00 H \nATOM 16 HG23 ILE A 1 -10.867 5.721 -5.364 1.00 0.00 H \nATOM 17 HD11 ILE A 1 -14.329 6.679 -8.366 1.00 0.00 H \nATOM 18 HD12 ILE A 1 -12.740 6.030 -7.962 1.00 0.00 H \nATOM 19 HD13 ILE A 1 -13.101 7.731 -7.663 1.00 0.00 H \nATOM 20 N ARG A 2 -12.031 5.020 -2.238 1.00 0.00 N \nATOM 21 CA ARG A 2 -10.996 4.938 -1.228 1.00 0.00 C \nATOM 22 C ARG A 2 -10.553 3.495 -1.030 1.00 0.00 C \nATOM 23 O ARG A 2 -11.381 2.588 -0.960 1.00 0.00 O \nATOM 24 CB ARG A 2 -11.490 5.514 0.107 1.00 0.00 C \nATOM 25 CG ARG A 2 -10.465 5.427 1.230 1.00 0.00 C \nATOM 26 CD ARG A 2 -9.262 6.317 0.963 1.00 0.00 C \nATOM 27 NE ARG A 2 -8.118 5.956 1.797 1.00 0.00 N \nATOM 28 CZ ARG A 2 -6.983 6.648 1.864 1.00 0.00 C \nATOM 29 NH1 ARG A 2 -6.846 7.782 1.184 1.00 0.00 N \nATOM 30 NH2 ARG A 2 -5.987 6.202 2.615 1.00 0.00 N \nATOM 31 H ARG A 2 -12.457 4.195 -2.540 1.00 0.00 H \nATOM 32 HA ARG A 2 -10.157 5.515 -1.572 1.00 0.00 H \nATOM 33 HB2 ARG A 2 -11.747 6.553 -0.036 1.00 0.00 H \nATOM 34 HB3 ARG A 2 -12.374 4.974 0.414 1.00 0.00 H \nATOM 35 HG2 ARG A 2 -10.930 5.737 2.153 1.00 0.00 H \nATOM 36 HG3 ARG A 2 -10.131 4.403 1.319 1.00 0.00 H \nATOM 37 HD2 ARG A 2 -8.982 6.220 -0.077 1.00 0.00 H \nATOM 38 HD3 ARG A 2 -9.536 7.342 1.166 1.00 0.00 H \nATOM 39 HE ARG A 2 -8.191 5.118 2.334 1.00 0.00 H \nATOM 40 HH11 ARG A 2 -7.596 8.121 0.617 1.00 0.00 H \nATOM 41 HH12 ARG A 2 -5.990 8.299 1.235 1.00 0.00 H \nATOM 42 HH21 ARG A 2 -6.088 5.334 3.133 1.00 0.00 H \nATOM 43 HH22 ARG A 2 -5.127 6.706 2.677 1.00 0.00 H \nATOM 44 N ASP A 3 -9.246 3.307 -0.925 1.00 0.00 N \nATOM 45 CA ASP A 3 -8.653 1.995 -0.711 1.00 0.00 C \nATOM 46 C ASP A 3 -7.319 2.165 -0.007 1.00 0.00 C \nATOM 47 O ASP A 3 -6.489 2.980 -0.425 1.00 0.00 O \nATOM 48 CB ASP A 3 -8.442 1.253 -2.032 1.00 0.00 C \nATOM 49 CG ASP A 3 -7.601 0.001 -1.854 1.00 0.00 C \nATOM 50 OD1 ASP A 3 -8.029 -0.906 -1.107 1.00 0.00 O \nATOM 51 OD2 ASP A 3 -6.502 -0.065 -2.443 1.00 0.00 O \nATOM 52 H ASP A 3 -8.652 4.081 -0.978 1.00 0.00 H \nATOM 53 HA ASP A 3 -9.317 1.423 -0.080 1.00 0.00 H \nATOM 54 HB2 ASP A 3 -9.401 0.969 -2.438 1.00 0.00 H \nATOM 55 HB3 ASP A 3 -7.939 1.908 -2.731 1.00 0.00 H \nHETATM 56 N CGU A 4 -7.117 1.414 1.061 1.00 0.00 N \nHETATM 57 CA CGU A 4 -5.883 1.503 1.816 1.00 0.00 C \nHETATM 58 C CGU A 4 -4.825 0.568 1.248 1.00 0.00 C \nHETATM 59 O CGU A 4 -4.335 -0.330 1.933 1.00 0.00 O \nHETATM 60 CB CGU A 4 -6.123 1.219 3.298 1.00 0.00 C \nHETATM 61 CG CGU A 4 -5.865 2.430 4.184 1.00 0.00 C \nHETATM 62 CD1 CGU A 4 -7.092 3.333 4.198 1.00 0.00 C \nHETATM 63 CD2 CGU A 4 -4.516 3.045 3.826 1.00 0.00 C \nHETATM 64 OE11 CGU A 4 -7.550 3.743 3.108 1.00 0.00 O \nHETATM 65 OE12 CGU A 4 -7.605 3.624 5.295 1.00 0.00 O \nHETATM 66 OE21 CGU A 4 -4.446 4.273 3.622 1.00 0.00 O \nHETATM 67 OE22 CGU A 4 -3.525 2.287 3.734 1.00 0.00 O \nHETATM 68 H CGU A 4 -7.815 0.792 1.350 1.00 0.00 H \nHETATM 69 HA CGU A 4 -5.523 2.518 1.717 1.00 0.00 H \nHETATM 70 HB2 CGU A 4 -5.467 0.423 3.612 1.00 0.00 H \nHETATM 71 HB3 CGU A 4 -7.149 0.911 3.435 1.00 0.00 H \nHETATM 72 HG CGU A 4 -5.762 2.055 5.191 1.00 0.00 H \nATOM 73 N CYS A 5 -4.475 0.798 -0.010 1.00 0.00 N \nATOM 74 CA CYS A 5 -3.472 0.008 -0.700 1.00 0.00 C \nATOM 75 C CYS A 5 -2.121 0.196 -0.002 1.00 0.00 C \nATOM 76 O CYS A 5 -1.236 -0.649 -0.079 1.00 0.00 O \nATOM 77 CB CYS A 5 -3.408 0.469 -2.164 1.00 0.00 C \nATOM 78 SG CYS A 5 -2.963 -0.811 -3.398 1.00 0.00 S \nATOM 79 H CYS A 5 -4.906 1.534 -0.497 1.00 0.00 H \nATOM 80 HA CYS A 5 -3.762 -1.032 -0.658 1.00 0.00 H \nATOM 81 HB2 CYS A 5 -4.376 0.853 -2.443 1.00 0.00 H \nATOM 82 HB3 CYS A 5 -2.686 1.266 -2.238 1.00 0.00 H \nATOM 83 N CYS A 6 -1.986 1.317 0.702 1.00 0.00 N \nATOM 84 CA CYS A 6 -0.766 1.632 1.431 1.00 0.00 C \nATOM 85 C CYS A 6 -0.529 0.619 2.552 1.00 0.00 C \nATOM 86 O CYS A 6 0.538 0.009 2.638 1.00 0.00 O \nATOM 87 CB CYS A 6 -0.867 3.047 2.007 1.00 0.00 C \nATOM 88 SG CYS A 6 0.734 3.813 2.414 1.00 0.00 S \nATOM 89 H CYS A 6 -2.737 1.944 0.743 1.00 0.00 H \nATOM 90 HA CYS A 6 0.060 1.590 0.738 1.00 0.00 H \nATOM 91 HB2 CYS A 6 -1.360 3.683 1.287 1.00 0.00 H \nATOM 92 HB3 CYS A 6 -1.458 3.015 2.910 1.00 0.00 H \nATOM 93 N SER A 7 -1.534 0.446 3.403 1.00 0.00 N \nATOM 94 CA SER A 7 -1.449 -0.488 4.518 1.00 0.00 C \nATOM 95 C SER A 7 -1.678 -1.926 4.050 1.00 0.00 C \nATOM 96 O SER A 7 -1.250 -2.879 4.706 1.00 0.00 O \nATOM 97 CB SER A 7 -2.476 -0.114 5.588 1.00 0.00 C \nATOM 98 OG SER A 7 -2.348 1.248 5.962 1.00 0.00 O \nATOM 99 H SER A 7 -2.362 0.968 3.281 1.00 0.00 H \nATOM 100 HA SER A 7 -0.459 -0.412 4.941 1.00 0.00 H \nATOM 101 HB2 SER A 7 -3.472 -0.276 5.202 1.00 0.00 H \nATOM 102 HB3 SER A 7 -2.325 -0.730 6.462 1.00 0.00 H \nATOM 103 HG SER A 7 -2.794 1.804 5.302 1.00 0.00 H \nATOM 104 N ASN A 8 -2.372 -2.066 2.922 1.00 0.00 N \nATOM 105 CA ASN A 8 -2.693 -3.372 2.346 1.00 0.00 C \nATOM 106 C ASN A 8 -1.437 -4.227 2.182 1.00 0.00 C \nATOM 107 O ASN A 8 -0.510 -3.851 1.468 1.00 0.00 O \nATOM 108 CB ASN A 8 -3.359 -3.187 0.980 1.00 0.00 C \nATOM 109 CG ASN A 8 -4.042 -4.442 0.468 1.00 0.00 C \nATOM 110 OD1 ASN A 8 -3.427 -5.503 0.360 1.00 0.00 O \nATOM 111 ND2 ASN A 8 -5.316 -4.321 0.125 1.00 0.00 N \nATOM 112 H ASN A 8 -2.692 -1.256 2.465 1.00 0.00 H \nATOM 113 HA ASN A 8 -3.379 -3.874 3.010 1.00 0.00 H \nATOM 114 HB2 ASN A 8 -4.095 -2.405 1.052 1.00 0.00 H \nATOM 115 HB3 ASN A 8 -2.607 -2.894 0.262 1.00 0.00 H \nATOM 116 HD21 ASN A 8 -5.742 -3.438 0.218 1.00 0.00 H \nATOM 117 HD22 ASN A 8 -5.781 -5.112 -0.213 1.00 0.00 H \nATOM 118 N PRO A 9 -1.396 -5.394 2.845 1.00 0.00 N \nATOM 119 CA PRO A 9 -0.250 -6.302 2.773 1.00 0.00 C \nATOM 120 C PRO A 9 0.021 -6.782 1.351 1.00 0.00 C \nATOM 121 O PRO A 9 1.165 -6.774 0.895 1.00 0.00 O \nATOM 122 CB PRO A 9 -0.650 -7.476 3.675 1.00 0.00 C \nATOM 123 CG PRO A 9 -2.130 -7.374 3.824 1.00 0.00 C \nATOM 124 CD PRO A 9 -2.459 -5.914 3.719 1.00 0.00 C \nATOM 125 HA PRO A 9 0.643 -5.833 3.163 1.00 0.00 H \nATOM 126 HB2 PRO A 9 -0.365 -8.406 3.204 1.00 0.00 H \nATOM 127 HB3 PRO A 9 -0.152 -7.385 4.629 1.00 0.00 H \nATOM 128 HG2 PRO A 9 -2.618 -7.926 3.034 1.00 0.00 H \nATOM 129 HG3 PRO A 9 -2.429 -7.757 4.789 1.00 0.00 H \nATOM 130 HD2 PRO A 9 -3.432 -5.777 3.268 1.00 0.00 H \nATOM 131 HD3 PRO A 9 -2.422 -5.447 4.691 1.00 0.00 H \nATOM 132 N ALA A 10 -1.034 -7.186 0.651 1.00 0.00 N \nATOM 133 CA ALA A 10 -0.908 -7.659 -0.720 1.00 0.00 C \nATOM 134 C ALA A 10 -0.405 -6.541 -1.622 1.00 0.00 C \nATOM 135 O ALA A 10 0.514 -6.738 -2.416 1.00 0.00 O \nATOM 136 CB ALA A 10 -2.239 -8.195 -1.226 1.00 0.00 C \nATOM 137 H ALA A 10 -1.922 -7.154 1.065 1.00 0.00 H \nATOM 138 HA ALA A 10 -0.190 -8.467 -0.730 1.00 0.00 H \nATOM 139 HB1 ALA A 10 -2.188 -8.333 -2.296 1.00 0.00 H \nATOM 140 HB2 ALA A 10 -3.024 -7.492 -0.990 1.00 0.00 H \nATOM 141 HB3 ALA A 10 -2.449 -9.142 -0.751 1.00 0.00 H \nATOM 142 N CYS A 11 -0.999 -5.362 -1.472 1.00 0.00 N \nATOM 143 CA CYS A 11 -0.604 -4.199 -2.255 1.00 0.00 C \nATOM 144 C CYS A 11 0.855 -3.860 -1.985 1.00 0.00 C \nATOM 145 O CYS A 11 1.616 -3.588 -2.908 1.00 0.00 O \nATOM 146 CB CYS A 11 -1.497 -3.004 -1.924 1.00 0.00 C \nATOM 147 SG CYS A 11 -1.125 -1.492 -2.874 1.00 0.00 S \nATOM 148 H CYS A 11 -1.715 -5.270 -0.806 1.00 0.00 H \nATOM 149 HA CYS A 11 -0.717 -4.447 -3.299 1.00 0.00 H \nATOM 150 HB2 CYS A 11 -2.524 -3.269 -2.123 1.00 0.00 H \nATOM 151 HB3 CYS A 11 -1.390 -2.767 -0.876 1.00 0.00 H \nATOM 152 N ARG A 12 1.236 -3.899 -0.715 1.00 0.00 N \nATOM 153 CA ARG A 12 2.610 -3.615 -0.305 1.00 0.00 C \nATOM 154 C ARG A 12 3.584 -4.605 -0.944 1.00 0.00 C \nATOM 155 O ARG A 12 4.677 -4.226 -1.357 1.00 0.00 O \nATOM 156 CB ARG A 12 2.727 -3.676 1.219 1.00 0.00 C \nATOM 157 CG ARG A 12 4.119 -3.357 1.738 1.00 0.00 C \nATOM 158 CD ARG A 12 4.195 -3.482 3.253 1.00 0.00 C \nATOM 159 NE ARG A 12 3.794 -4.812 3.732 1.00 0.00 N \nATOM 160 CZ ARG A 12 4.448 -5.950 3.464 1.00 0.00 C \nATOM 161 NH1 ARG A 12 5.551 -5.936 2.720 1.00 0.00 N \nATOM 162 NH2 ARG A 12 3.994 -7.103 3.947 1.00 0.00 N \nATOM 163 H ARG A 12 0.570 -4.133 -0.026 1.00 0.00 H \nATOM 164 HA ARG A 12 2.864 -2.618 -0.641 1.00 0.00 H \nATOM 165 HB2 ARG A 12 2.035 -2.969 1.652 1.00 0.00 H \nATOM 166 HB3 ARG A 12 2.464 -4.670 1.549 1.00 0.00 H \nATOM 167 HG2 ARG A 12 4.823 -4.046 1.295 1.00 0.00 H \nATOM 168 HG3 ARG A 12 4.376 -2.347 1.456 1.00 0.00 H \nATOM 169 HD2 ARG A 12 5.212 -3.294 3.564 1.00 0.00 H \nATOM 170 HD3 ARG A 12 3.545 -2.741 3.694 1.00 0.00 H \nATOM 171 HE ARG A 12 2.983 -4.853 4.287 1.00 0.00 H \nATOM 172 HH11 ARG A 12 5.902 -5.070 2.355 1.00 0.00 H \nATOM 173 HH12 ARG A 12 6.036 -6.790 2.516 1.00 0.00 H \nATOM 174 HH21 ARG A 12 3.166 -7.122 4.512 1.00 0.00 H \nATOM 175 HH22 ARG A 12 4.478 -7.958 3.751 1.00 0.00 H \nATOM 176 N VAL A 13 3.185 -5.870 -1.025 1.00 0.00 N \nATOM 177 CA VAL A 13 4.034 -6.901 -1.623 1.00 0.00 C \nATOM 178 C VAL A 13 4.150 -6.701 -3.134 1.00 0.00 C \nATOM 179 O VAL A 13 5.225 -6.847 -3.711 1.00 0.00 O \nATOM 180 CB VAL A 13 3.500 -8.322 -1.325 1.00 0.00 C \nATOM 181 CG1 VAL A 13 4.302 -9.378 -2.073 1.00 0.00 C \nATOM 182 CG2 VAL A 13 3.534 -8.600 0.170 1.00 0.00 C \nATOM 183 H VAL A 13 2.300 -6.118 -0.677 1.00 0.00 H \nATOM 184 HA VAL A 13 5.017 -6.811 -1.187 1.00 0.00 H \nATOM 185 HB VAL A 13 2.474 -8.378 -1.656 1.00 0.00 H \nATOM 186 HG11 VAL A 13 4.872 -9.965 -1.368 1.00 0.00 H \nATOM 187 HG12 VAL A 13 4.976 -8.894 -2.765 1.00 0.00 H \nATOM 188 HG13 VAL A 13 3.629 -10.022 -2.618 1.00 0.00 H \nATOM 189 HG21 VAL A 13 3.140 -9.587 0.362 1.00 0.00 H \nATOM 190 HG22 VAL A 13 2.933 -7.866 0.687 1.00 0.00 H \nATOM 191 HG23 VAL A 13 4.552 -8.544 0.524 1.00 0.00 H \nATOM 192 N ASN A 14 3.035 -6.367 -3.768 1.00 0.00 N \nATOM 193 CA ASN A 14 3.003 -6.148 -5.210 1.00 0.00 C \nATOM 194 C ASN A 14 3.699 -4.849 -5.578 1.00 0.00 C \nATOM 195 O ASN A 14 4.495 -4.800 -6.514 1.00 0.00 O \nATOM 196 CB ASN A 14 1.562 -6.086 -5.702 1.00 0.00 C \nATOM 197 CG ASN A 14 0.807 -7.397 -5.550 1.00 0.00 C \nATOM 198 OD1 ASN A 14 -0.395 -7.462 -5.800 1.00 0.00 O \nATOM 199 ND2 ASN A 14 1.500 -8.454 -5.148 1.00 0.00 N \nATOM 200 H ASN A 14 2.206 -6.266 -3.253 1.00 0.00 H \nATOM 201 HA ASN A 14 3.508 -6.970 -5.689 1.00 0.00 H \nATOM 202 HB2 ASN A 14 1.047 -5.323 -5.145 1.00 0.00 H \nATOM 203 HB3 ASN A 14 1.564 -5.815 -6.745 1.00 0.00 H \nATOM 204 HD21 ASN A 14 2.457 -8.340 -4.969 1.00 0.00 H \nATOM 205 HD22 ASN A 14 1.026 -9.304 -5.046 1.00 0.00 H \nATOM 206 N ASN A 15 3.382 -3.798 -4.836 1.00 0.00 N \nATOM 207 CA ASN A 15 3.960 -2.480 -5.073 1.00 0.00 C \nATOM 208 C ASN A 15 4.714 -1.982 -3.840 1.00 0.00 C \nATOM 209 O ASN A 15 4.149 -1.291 -2.989 1.00 0.00 O \nATOM 210 CB ASN A 15 2.857 -1.479 -5.439 1.00 0.00 C \nATOM 211 CG ASN A 15 2.056 -1.890 -6.659 1.00 0.00 C \nATOM 212 OD1 ASN A 15 2.587 -2.000 -7.761 1.00 0.00 O \nATOM 213 ND2 ASN A 15 0.763 -2.117 -6.467 1.00 0.00 N \nATOM 214 H ASN A 15 2.734 -3.908 -4.105 1.00 0.00 H \nATOM 215 HA ASN A 15 4.653 -2.566 -5.894 1.00 0.00 H \nATOM 216 HB2 ASN A 15 2.177 -1.386 -4.605 1.00 0.00 H \nATOM 217 HB3 ASN A 15 3.308 -0.517 -5.635 1.00 0.00 H \nATOM 218 HD21 ASN A 15 0.402 -2.007 -5.561 1.00 0.00 H \nATOM 219 HD22 ASN A 15 0.224 -2.386 -7.238 1.00 0.00 H \nHETATM 220 N HYP A 16 6.007 -2.327 -3.718 1.00 0.00 N \nHETATM 221 CA HYP A 16 6.820 -1.910 -2.578 1.00 0.00 C \nHETATM 222 C HYP A 16 7.263 -0.454 -2.655 1.00 0.00 C \nHETATM 223 O HYP A 16 7.793 -0.007 -3.670 1.00 0.00 O \nHETATM 224 CB HYP A 16 8.031 -2.830 -2.666 1.00 0.00 C \nHETATM 225 CG HYP A 16 7.504 -4.080 -3.281 1.00 0.00 C \nHETATM 226 CD2 HYP A 16 6.297 -3.718 -4.107 1.00 0.00 C \nHETATM 227 OD1 HYP A 16 7.105 -5.002 -2.289 1.00 0.00 O \nHETATM 228 HA HYP A 16 6.295 -2.080 -1.645 1.00 0.00 H \nHETATM 229 HB1 HYP A 16 8.424 -3.011 -1.677 1.00 0.00 H \nHETATM 230 HB2 HYP A 16 8.790 -2.372 -3.285 1.00 0.00 H \nHETATM 231 HG HYP A 16 8.274 -4.579 -3.866 1.00 0.00 H \nHETATM 232 HD21 HYP A 16 6.516 -3.786 -5.161 1.00 0.00 H \nHETATM 233 HD22 HYP A 16 5.464 -4.359 -3.865 1.00 0.00 H \nHETATM 234 HD1 HYP A 16 6.161 -4.854 -2.153 1.00 0.00 H \nATOM 235 N HIS A 17 7.059 0.260 -1.552 1.00 0.00 N \nATOM 236 CA HIS A 17 7.438 1.666 -1.422 1.00 0.00 C \nATOM 237 C HIS A 17 6.795 2.530 -2.500 1.00 0.00 C \nATOM 238 O HIS A 17 7.308 3.586 -2.859 1.00 0.00 O \nATOM 239 CB HIS A 17 8.966 1.803 -1.441 1.00 0.00 C \nATOM 240 CG HIS A 17 9.485 3.144 -0.996 1.00 0.00 C \nATOM 241 ND1 HIS A 17 10.804 3.524 -1.140 1.00 0.00 N \nATOM 242 CD2 HIS A 17 8.861 4.188 -0.392 1.00 0.00 C \nATOM 243 CE1 HIS A 17 10.968 4.737 -0.645 1.00 0.00 C \nATOM 244 NE2 HIS A 17 9.807 5.162 -0.186 1.00 0.00 N \nATOM 245 H HIS A 17 6.652 -0.180 -0.785 1.00 0.00 H \nATOM 246 HA HIS A 17 7.068 2.009 -0.472 1.00 0.00 H \nATOM 247 HB2 HIS A 17 9.391 1.056 -0.787 1.00 0.00 H \nATOM 248 HB3 HIS A 17 9.313 1.626 -2.445 1.00 0.00 H \nATOM 249 HD1 HIS A 17 11.512 2.986 -1.550 1.00 0.00 H \nATOM 250 HD2 HIS A 17 7.815 4.241 -0.123 1.00 0.00 H \nATOM 251 HE1 HIS A 17 11.896 5.290 -0.621 1.00 0.00 H \nATOM 252 N VAL A 18 5.647 2.098 -2.976 1.00 0.00 N \nATOM 253 CA VAL A 18 4.915 2.858 -3.962 1.00 0.00 C \nATOM 254 C VAL A 18 4.167 3.981 -3.251 1.00 0.00 C \nATOM 255 O VAL A 18 4.044 5.095 -3.755 1.00 0.00 O \nATOM 256 CB VAL A 18 3.949 1.957 -4.757 1.00 0.00 C \nATOM 257 CG1 VAL A 18 2.729 2.731 -5.214 1.00 0.00 C \nATOM 258 CG2 VAL A 18 4.671 1.349 -5.951 1.00 0.00 C \nATOM 259 H VAL A 18 5.267 1.263 -2.632 1.00 0.00 H \nATOM 260 HA VAL A 18 5.629 3.292 -4.648 1.00 0.00 H \nATOM 261 HB VAL A 18 3.622 1.153 -4.115 1.00 0.00 H \nATOM 262 HG11 VAL A 18 2.422 3.398 -4.422 1.00 0.00 H \nATOM 263 HG12 VAL A 18 1.929 2.044 -5.440 1.00 0.00 H \nATOM 264 HG13 VAL A 18 2.975 3.305 -6.093 1.00 0.00 H \nATOM 265 HG21 VAL A 18 3.967 1.187 -6.753 1.00 0.00 H \nATOM 266 HG22 VAL A 18 5.112 0.407 -5.663 1.00 0.00 H \nATOM 267 HG23 VAL A 18 5.447 2.023 -6.283 1.00 0.00 H \nATOM 268 N CYS A 19 3.702 3.671 -2.051 1.00 0.00 N \nATOM 269 CA CYS A 19 2.993 4.639 -1.228 1.00 0.00 C \nATOM 270 C CYS A 19 3.979 5.407 -0.355 1.00 0.00 C \nATOM 271 O CYS A 19 5.079 4.918 -0.073 1.00 0.00 O \nATOM 272 CB CYS A 19 1.960 3.929 -0.354 1.00 0.00 C \nATOM 273 SG CYS A 19 1.099 5.021 0.823 1.00 0.00 S \nATOM 274 H CYS A 19 3.861 2.772 -1.698 1.00 0.00 H \nATOM 275 HA CYS A 19 2.489 5.332 -1.884 1.00 0.00 H \nATOM 276 HB2 CYS A 19 1.212 3.479 -0.987 1.00 0.00 H \nATOM 277 HB3 CYS A 19 2.452 3.156 0.217 1.00 0.00 H \nTER 278 CYS A 19 \n';Mol1.loadMolecule(ChemDoodle.readPDB(Fmol, 1));Mol1.startAnimation();Mol1.stopAnimation();var resType = 'amino';function setDisp1(obj){Mol1.specs["proteins_display"+obj.value]=obj.checked;Mol1.repaint()}function setProj1(yesPers){Mol1.specs.projectionPerspective_3D=yesPers;Mol1.setupScene();Mol1.repaint()}function setColor1(obj){Mol1.specs.proteins_residueColor=resType=obj.value;Mol1.repaint()}function resColor1(obj){this.checked ? Mol1.specs.proteins_residueColor=resType : Mol1.specs.proteins_residueColor='none';Mol1.repaint()}</script> 视图: <input type="radio" name="group2" onclick="setProj1(true)">投影 <input type="radio" name="group2" onclick="setProj1(false)" checked="">正交 着色: <input type="checkbox" onclick="Mol1.specs.macro_colorByChain=this.checked;Mol1.repaint()">按链 <input type="checkbox" onclick="resColor1(this)">按残基 模式: <input type="checkbox" value="Ribbon" checked="" onclick="setDisp1(this)">飘带 <input type="checkbox" value="Backbone" onclick="setDisp1(this)">骨架 <input type="checkbox" value="PipePlank" onclick="setDisp1(this)">管板 <input type="checkbox" checked="" onclick="Mol1.specs.proteins_ribbonCartoonize=this.checked;Mol1.repaint()">卡通 显示: <input type="checkbox" onclick="Mol1.specs.macro_showWater=this.checked;Mol1.repaint()">水分子 <input type="checkbox" onclick="Mol1.specs.atoms_nonBondedAsStars_3D=this.checked;Mol1.repaint()">非键原子 <input type="checkbox" onclick="Mol1.specs.atoms_displayLabels_3D=this.checked;Mol1.repaint()">名称 颜色: <input type="radio" value="amino" name="radio2" onclick="setColor1(this)">氨基酸 <input type="radio" value="shapely" name="radio2" onclick="setColor1(this)">形状 <input type="radio" value="polarity" name="radio2" onclick="setColor1(this)">极性 <input type="radio" value="acidity" name="radio2" onclick="setColor1(this)">酸性 <input type="radio" value="rainbow" name="radio2" onclick="setColor1(this)">彩虹 左键: 转动滚轮: 缩放双击: 自动旋转开关 Alt+左键: 移动 <figurecaption>Fig.1</figurecaption></figure>

查看pdb文件或通过可视化软件显示可知, GID中实际包含了两种非标准残基, CGU和HYP. 但GROMACS自带的amber99sb-ildn力场中包含了HYP残基的数据, 所以我们不需要对HYP使用非标准残基的处理流程, 只要稍加处理即可(具体处理见后文). 而对于CGU残基, 我们必须按照非标准残基的流程进行处理. 幸运的是, HYP残基为我们提供了一个非标准残基的示例, 如果在处理CGU时遇到问题(主要是准备其rtp文件), 可以参考amber99sb-ildn力场中HYP的数据.

第二步更新残基类型文件`residuetypes.dat`

简单地讲, 对GROMACS而言, 所谓非标准残基就是指pdb2gmx不能识别的化学基团, 或者说是GROMACS的残基数据库中没有对应数据的化学基团. 为了使pdb2gmx能够识别新定义的残基, 并进而正确处理, 我们首先要将新残基的名称添加到/GMX安装路径/share/gromacs/top/residuetypes.dat文件中, 并指定新残基的类型. 浏览下此文件就可以看出, 残基的可能类型有Protein(蛋白), DNA, RNA, Water(水), Ion(离子). 我们要添加的CGU残基属于蛋白, 且处于肽链中间, 而不处于肽链两端, 所以只要简单地在residuetypes.dat文件中添加一行指明其类型是蛋白就可以了. 为简单起见, 我们将其添加到文件的最前面, 如下所示

<table class="highlighttable"><th colspan="2" style="text-align:left">residuetypes.dat</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">CGU Protein ; 新添加的CGU残基 ABU Protein ACE Protein AIB Protein ALA Protein ARG Protein ARGN Protein ASN Protein ASN1 Protein ASP Protein ASP1 Protein ...(下略) </pre></div> </td></tr></table>

这里多说几句, 如果你浏览一下residuetypes.dat文件的内容, 可能会发现

如果要添加的蛋白残基有可能处于端基位置, 还需要定义它处于N端和/或C端时的残基名称
如果要添加的残基属于DNA或RNA, 且有可能处于3端和/或5端, 同样需要定义其处于端基时的名称
如果要添加的残基可能存在不同的质子化状态, 同样需要定义其处于不同质子化状态时的名称

具体如何定义, 相信你能够从文件中得到答案.

第三步创建包含残基的完整分子

为了将CGU残基的拓扑信息添加到GROMACS残基数据库中, 我们首先需要得到CGU的拓扑数据. 为此, 我们需要先将CGU残基从肽链中截取下来, 形成一个包含CGU的, 单独的, 完整的分子, 从而能够利用工具获得它的拓扑文件. 由于我们使用的是AMBER力场, 所以我们将使用ambertools获得这个分子的拓扑文件.

在抽取分子时, 如果我们只截取属于CGU的原子, 那么得到的并不是一个完整的分子, 需要进行封端处理. 常规的封端方法是使用氢原子. 但对于我们的情况并不合适, 因为如果只是对CGU用氢原子进行封端形成完整分子, 那么ambertools在判断CGU中每个原子的类型时, 很可能会与实际情况不符, 因为在实际的肽链中, 与CGU前后相连的是氨基酸残基的碳原子和氮原子. 所以最好是将CGU前后相邻的两个残基一起截取下来, 然后将两个相邻的残基用氢原子封端. 这样做时, 如果得到的分子过大, 导致ambertools不易处理, 也可以根据情况, 只截取相邻残基中与CGU的氮原子和碳原子最接近的一部分原子, 再使用氢原子封端. 有时候, 也可以分别使用乙酰基和甲胺基对截取残基的氮原子和碳原子进行封端. 这应该是能够使用的最小的封端基团了.

我们在这里将截取与CGU相邻的两个残基, 3ASP与5CYS, 然后用氢原子对它们封端. 所以, 我们截取gid.pdb中属于3ASP-4CGU-5CYS残基的原子. 截取原子时, 我们可以手动操作, 也可以使用PyMOL等可视化软件. 将截取的原子坐标保存为cgu.pdb. 其内容如下

<table class="highlighttable"><th colspan="2" style="text-align:left">cgu.pdb</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">ATOM 44 N ASP A 3 -9.246 3.307 -0.925 1.00 0.00 N ATOM 45 CA ASP A 3 -8.653 1.995 -0.711 1.00 0.00 C ATOM 46 C ASP A 3 -7.319 2.165 -0.007 1.00 0.00 C ATOM 47 O ASP A 3 -6.489 2.980 -0.425 1.00 0.00 O ATOM 48 CB ASP A 3 -8.442 1.253 -2.032 1.00 0.00 C ATOM 49 CG ASP A 3 -7.601 0.001 -1.854 1.00 0.00 C ATOM 50 OD1 ASP A 3 -8.029 -0.906 -1.107 1.00 0.00 O ATOM 51 OD2 ASP A 3 -6.502 -0.065 -2.443 1.00 0.00 O ATOM 52 H ASP A 3 -8.652 4.081 -0.978 1.00 0.00 H ATOM 53 HA ASP A 3 -9.317 1.423 -0.080 1.00 0.00 H ATOM 54 HB2 ASP A 3 -9.401 0.969 -2.438 1.00 0.00 H ATOM 55 HB3 ASP A 3 -7.939 1.908 -2.731 1.00 0.00 H HETATM 56 N CGU A 4 -7.117 1.414 1.061 1.00 0.00 N HETATM 57 CA CGU A 4 -5.883 1.503 1.816 1.00 0.00 C HETATM 58 C CGU A 4 -4.825 0.568 1.248 1.00 0.00 C HETATM 59 O CGU A 4 -4.335 -0.330 1.933 1.00 0.00 O HETATM 60 CB CGU A 4 -6.123 1.219 3.298 1.00 0.00 C HETATM 61 CG CGU A 4 -5.865 2.430 4.184 1.00 0.00 C HETATM 62 CD1 CGU A 4 -7.092 3.333 4.198 1.00 0.00 C HETATM 63 CD2 CGU A 4 -4.516 3.045 3.826 1.00 0.00 C HETATM 64 OE11 CGU A 4 -7.550 3.743 3.108 1.00 0.00 O HETATM 65 OE12 CGU A 4 -7.605 3.624 5.295 1.00 0.00 O HETATM 66 OE21 CGU A 4 -4.446 4.273 3.622 1.00 0.00 O HETATM 67 OE22 CGU A 4 -3.525 2.287 3.734 1.00 0.00 O HETATM 68 H CGU A 4 -7.815 0.792 1.350 1.00 0.00 H HETATM 69 HA CGU A 4 -5.523 2.518 1.717 1.00 0.00 H HETATM 70 HB2 CGU A 4 -5.467 0.423 3.612 1.00 0.00 H HETATM 71 HB3 CGU A 4 -7.149 0.911 3.435 1.00 0.00 H HETATM 72 HG CGU A 4 -5.762 2.055 5.191 1.00 0.00 H ATOM 73 N CYS A 5 -4.475 0.798 -0.010 1.00 0.00 N ATOM 74 CA CYS A 5 -3.472 0.008 -0.700 1.00 0.00 C ATOM 75 C CYS A 5 -2.121 0.196 -0.002 1.00 0.00 C ATOM 76 O CYS A 5 -1.236 -0.649 -0.079 1.00 0.00 O ATOM 77 CB CYS A 5 -3.408 0.469 -2.164 1.00 0.00 C ATOM 78 SG CYS A 5 -2.963 -0.811 -3.398 1.00 0.00 S ATOM 79 H CYS A 5 -4.906 1.534 -0.497 1.00 0.00 H ATOM 80 HA CYS A 5 -3.762 -1.032 -0.658 1.00 0.00 H ATOM 81 HB2 CYS A 5 -4.376 0.853 -2.443 1.00 0.00 H ATOM 82 HB3 CYS A 5 -2.686 1.266 -2.238 1.00 0.00 H </pre></div> </td></tr></table>

使用分子编辑软件(我使用的是GaussView)打开cgu.pdb, 检查下分子, 看是否需要补加缺失的原子(大多是氢原子).

<figure><script>var Mol2=new ChemDoodle.TransformCanvas3D('Mol-2',642,396);Mol2.specs.shapes_color='#fff';Mol2.specs.backgroundColor='black';Mol2.specs.set3DRepresentation('Ball and Stick');Mol2.specs.projectionPerspective_3D=false;Mol2.specs.compass_display=true; /*//Mol2.specs.atoms_resolution_3D=15; //Mol2.specs.bonds_resolution_3D=15; //Mol2.specs.crystals_unitCellLineWidth=1.5;*/ Mol2.nextFrame=function(delta){var matrix=[];ChemDoodle.lib.mat4.identity(matrix);var change=delta*Math.PI/15000;ChemDoodle.lib.mat4.rotate(matrix,change,[1,0,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,1,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,0,1]);ChemDoodle.lib.mat4.multiply(this.rotationMatrix, matrix)}; Mol2.startAnimation=ChemDoodle._AnimatorCanvas.prototype.startAnimation;Mol2.stopAnimation=ChemDoodle._AnimatorCanvas.prototype.stopAnimation;Mol2.isRunning=ChemDoodle._AnimatorCanvas.prototype.isRunning;Mol2.dblclick=ChemDoodle.RotatorCanvas.prototype.dblclick;Mol2.timeout=5;Mol2.handle=null; var Fmol='39\ncgu\nN -9.246 3.307 -0.925\nC -8.653 1.995 -0.711\nC -7.319 2.165 -0.007\nO -6.489 2.980 -0.425\nC -8.442 1.253 -2.032\nC -7.601 0.001 -1.854\nO -8.029 -0.906 -1.107\nO -6.502 -0.065 -2.443\nH -8.652 4.081 -0.978\nH -9.317 1.423 -0.080\nH -9.401 0.969 -2.438\nH -7.939 1.908 -2.731\nN -7.117 1.414 1.061\nC -5.883 1.503 1.816\nC -4.825 0.568 1.248\nO -4.335 -0.330 1.933\nC -6.123 1.219 3.298\nC -5.865 2.430 4.184\nC -7.092 3.333 4.198\nC -4.516 3.045 3.826\nO -7.550 3.743 3.108\nO -7.605 3.624 5.295\nO -4.446 4.273 3.622\nO -3.525 2.287 3.734\nH -7.815 0.792 1.350\nH -5.523 2.518 1.717\nH -5.467 0.423 3.612\nH -7.149 0.911 3.435\nH -5.762 2.055 5.191\nN -4.475 0.798 -0.010\nC -3.472 0.008 -0.700\nC -2.121 0.196 -0.002\nO -1.236 -0.649 -0.079\nC -3.408 0.469 -2.164\nS -2.963 -0.811 -3.398\nH -4.906 1.534 -0.497\nH -3.762 -1.032 -0.658\nH -4.376 0.853 -2.443\nH -2.686 1.266 -2.238\n'; Mol2.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol2.startAnimation();Mol2.stopAnimation();function setProj2(yesPers){Mol2.specs.projectionPerspective_3D=yesPers;Mol2.setupScene();Mol2.repaint()}function setModel2(model){Mol2.specs.set3DRepresentation(model);Mol2.setupScene();Mol2.repaint()}function setSpeed2(){Mol2.timeout=500-document.getElementById('spd2').value;Mol2.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol2.startAnimation()}</script> 视图: <input type='radio' name='group2' onclick='setProj2(true)'>投影 <input type='radio' name='group2' onclick='setProj2(false)' checked=''>正交速度: <input type='range' id='spd2' min='1' max='500' onchange='setSpeed2()'/> 模型: <input type='radio' name='model' onclick='setModel2('Ball and Stick')' checked=''>球棍 <input type='radio' name='model' onclick='setModel2('van der Waals Spheres')'>范德华球 <input type='radio' name='model' onclick='setModel2('Stick')'>棍状 <input type='radio' name='model' onclick='setModel2('Wireframe')'>线框 <input type='radio' name='model' onclick='setModel2('Line')'>线型 <input type='checkbox' onclick='Mol2.specs.atoms_displayLabels_3D=this.checked;Mol2.repaint()'>名称 左键: 转动滚轮: 缩放双击: 自动旋转开关 Alt+左键: 移动 <figurecaption>Fig.2</figurecaption></figure>

可以看到, CGU残基的两个羧基上都没有氢原子, 需要根据模拟条件(如pH)决定是否添加(也就是决定其质子化状态). 我们在这里假定CGU的O22和O24都需要补加氢原子. 此外, 3ASP的N1和O8, 5CYS的C32和S35也都需要使用氢原子进行封端. 使用GaussView的加氢功能补加这些氢原子, 然后将分子保存为cgu_h.pdb. GaussView给出的原子坐标如下

<table class="highlighttable"><th colspan="2" style="text-align:left">cgu_h.pdb</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">REMARK 1 File created by GaussView 5.0.9 ATOM 1 N ASP A 3 -9.246 3.307 -0.925 N ATOM 2 CA ASP A 3 -8.653 1.995 -0.711 C ATOM 3 C ASP A 3 -7.319 2.165 -0.007 C ATOM 4 O ASP A 3 -6.489 2.980 -0.425 O ATOM 5 CB ASP A 3 -8.442 1.253 -2.032 C ATOM 6 CG ASP A 3 -7.601 0.001 -1.854 C ATOM 7 OD1 ASP A 3 -8.029 -0.906 -1.107 O ATOM 8 OD2 ASP A 3 -6.502 -0.065 -2.443 O ATOM 9 H ASP A 3 -8.662 4.116 -0.983 H ATOM 10 HA ASP A 3 -9.317 1.423 -0.080 H ATOM 11 HB2 ASP A 3 -9.401 0.969 -2.438 H ATOM 12 HB3 ASP A 3 -7.939 1.908 -2.731 H HETATM 13 N CGU A 4 -7.117 1.414 1.061 N HETATM 14 CA CGU A 4 -5.883 1.503 1.816 C HETATM 15 C CGU A 4 -4.825 0.568 1.248 C HETATM 16 O CGU A 4 -4.335 -0.330 1.933 O HETATM 17 CB CGU A 4 -6.123 1.219 3.298 C HETATM 18 CG CGU A 4 -5.865 2.430 4.184 C HETATM 19 CD1 CGU A 4 -7.092 3.333 4.198 C HETATM 20 CD2 CGU A 4 -4.516 3.045 3.826 C HETATM 21 OE11 CGU A 4 -7.550 3.743 3.108 O HETATM 22 OE12 CGU A 4 -7.605 3.624 5.295 O HETATM 23 OE21 CGU A 4 -4.446 4.273 3.622 O HETATM 24 OE22 CGU A 4 -3.525 2.287 3.734 O HETATM 25 H CGU A 4 -7.815 0.792 1.350 H HETATM 26 HA CGU A 4 -5.523 2.518 1.717 H HETATM 27 HB2 CGU A 4 -5.467 0.423 3.612 H HETATM 28 HB3 CGU A 4 -7.149 0.911 3.435 H HETATM 29 HG CGU A 4 -5.762 2.055 5.191 H ATOM 30 N CYS A 5 -4.475 0.798 -0.010 N ATOM 31 CA CYS A 5 -3.472 0.008 -0.700 C ATOM 32 C CYS A 5 -2.121 0.196 -0.002 C ATOM 33 O CYS A 5 -1.233 -0.690 -0.101 O ATOM 34 CB CYS A 5 -3.408 0.469 -2.164 C ATOM 35 SG CYS A 5 -2.963 -0.811 -3.398 S ATOM 36 H CYS A 5 -4.906 1.534 -0.497 H ATOM 37 HA CYS A 5 -3.762 -1.032 -0.658 H ATOM 38 HB2 CYS A 5 -4.376 0.853 -2.443 H ATOM 39 HB3 CYS A 5 -2.686 1.266 -2.238 H TER 40 CYS A 5 HETATM 41 H 0 -6.912 3.794 5.937 H HETATM 42 H 0 -2.719 2.805 3.797 H HETATM 43 H 0 -10.238 3.399 -1.014 H HETATM 44 H 0 -5.791 0.007 -1.802 H HETATM 45 H 0 -1.932 1.081 0.569 H HETATM 46 H 0 -1.659 -0.900 -3.484 H </pre></div> </td></tr></table>

为了后面能够方便地使用ambertools获得这个分子的拓扑文件, 并整理拓扑信息, 我们需要对文件进行一些修改:

删除TER标识, 因为是单个分子
所有原子都使用HETATM类型, 因为不将分子视为肽链
补加的氢原子移动到其所属残基后面: 41, 42放到4CGU后面, 43, 44放到3ASP后面, 45, 46放到5CYS后面
调整残基位置, 4CGU放在最前面, 然后是3ASP, 最后是5CYS, 这样在后面整理拓扑信息时更方便
修改CGU中的氢原子的名称. 因为使用pdb2gmx自动加氢时, 对氢原子的名称有规定. 如果在一个原子上同时添加多个氢原子, 那么这些氢原子名称中的编号必须从1开始. 而CGU中HB2和HB3的编号不是从1开始, 所以要分别改为HB1和HB2. 同样, 补加的羧基氢原子分别更名为HE1和HE2.
将3ASP的氮原子和碳原子名称后面添加_, 5CYS的后面添加__. 这样可以区分不同残基, 方便后面整理拓扑信息.
所有原子的残基名称统一改为MOL, 链标识改为A, 残基编号为1.
3ASP和5CYS的氢原子名称可以不用管, 因为最终我们只需要CGU的拓扑信息.

注意:

GROMACS残基数据库中, 每个残基中各个原子的名称唯一且彼此间不能重复.
残基分子中与相邻氨基酸残基相连的原子名称必须为N和C, 否则形成拓扑文件时残基将无法与相邻残基形成正确的键.

将修改后的文件仍保存为cgu_H.pdb, 其内容如下:

<table class="highlighttable"><th colspan="2" style="text-align:left">cgu_H.pdb</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">HETATM 13 N MOL A 1 -7.117 1.414 1.061 N HETATM 14 CA MOL A 1 -5.883 1.503 1.816 C HETATM 15 C MOL A 1 -4.825 0.568 1.248 C HETATM 16 O MOL A 1 -4.335 -0.330 1.933 O HETATM 17 CB MOL A 1 -6.123 1.219 3.298 C HETATM 18 CG MOL A 1 -5.865 2.430 4.184 C HETATM 19 CD1 MOL A 1 -7.092 3.333 4.198 C HETATM 20 CD2 MOL A 1 -4.516 3.045 3.826 C HETATM 21 OE11 MOL A 1 -7.550 3.743 3.108 O HETATM 22 OE12 MOL A 1 -7.605 3.624 5.295 O HETATM 23 OE21 MOL A 1 -4.446 4.273 3.622 O HETATM 24 OE22 MOL A 1 -3.525 2.287 3.734 O HETATM 25 H MOL A 1 -7.815 0.792 1.350 H HETATM 26 HA MOL A 1 -5.523 2.518 1.717 H HETATM 27 HB1 MOL A 1 -5.467 0.423 3.612 H HETATM 28 HB2 MOL A 1 -7.149 0.911 3.435 H HETATM 29 HG MOL A 1 -5.762 2.055 5.191 H HETATM 41 HE1 MOL A 1 -6.912 3.794 5.937 H HETATM 42 HE2 MOL A 1 -2.719 2.805 3.797 H HETATM 1 N_ MOL A 1 -9.246 3.307 -0.925 N HETATM 2 CA MOL A 1 -8.653 1.995 -0.711 C HETATM 3 C_ MOL A 1 -7.319 2.165 -0.007 C HETATM 4 O MOL A 1 -6.489 2.980 -0.425 O HETATM 5 CB MOL A 1 -8.442 1.253 -2.032 C HETATM 6 CG MOL A 1 -7.601 0.001 -1.854 C HETATM 7 OD1 MOL A 1 -8.029 -0.906 -1.107 O HETATM 8 OD2 MOL A 1 -6.502 -0.065 -2.443 O HETATM 9 H MOL A 1 -8.662 4.116 -0.983 H HETATM 10 HA MOL A 1 -9.317 1.423 -0.080 H HETATM 11 HB2 MOL A 1 -9.401 0.969 -2.438 H HETATM 12 HB3 MOL A 1 -7.939 1.908 -2.731 H HETATM 43 H MOL A 1 -10.238 3.399 -1.014 H HETATM 44 H MOL A 1 -5.791 0.007 -1.802 H HETATM 30 N__ MOL A 1 -4.475 0.798 -0.010 N HETATM 31 CA MOL A 1 -3.472 0.008 -0.700 C HETATM 32 C__ MOL A 1 -2.121 0.196 -0.002 C HETATM 33 O MOL A 1 -1.233 -0.690 -0.101 O HETATM 34 CB MOL A 1 -3.408 0.469 -2.164 C HETATM 35 SG MOL A 1 -2.963 -0.811 -3.398 S HETATM 36 H MOL A 1 -4.906 1.534 -0.497 H HETATM 37 HA MOL A 1 -3.762 -1.032 -0.658 H HETATM 38 HB2 MOL A 1 -4.376 0.853 -2.443 H HETATM 39 HB3 MOL A 1 -2.686 1.266 -2.238 H HETATM 45 H MOL A 1 -1.932 1.081 0.569 H HETATM 46 H MOL A 1 -1.659 -0.900 -3.484 H </pre></div> </td></tr></table>

结构如下

<figure><script>var Mol3=new ChemDoodle.TransformCanvas3D('Mol-3',642,396);Mol3.specs.shapes_color='#fff';Mol3.specs.backgroundColor='black';Mol3.specs.set3DRepresentation('Ball and Stick');Mol3.specs.projectionPerspective_3D=false;Mol3.specs.compass_display=true; /*//Mol3.specs.atoms_resolution_3D=15; //Mol3.specs.bonds_resolution_3D=15; //Mol3.specs.crystals_unitCellLineWidth=1.5;*/ Mol3.nextFrame=function(delta){var matrix=[];ChemDoodle.lib.mat4.identity(matrix);var change=delta*Math.PI/15000;ChemDoodle.lib.mat4.rotate(matrix,change,[1,0,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,1,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,0,1]);ChemDoodle.lib.mat4.multiply(this.rotationMatrix, matrix)}; Mol3.startAnimation=ChemDoodle._AnimatorCanvas.prototype.startAnimation;Mol3.stopAnimation=ChemDoodle._AnimatorCanvas.prototype.stopAnimation;Mol3.isRunning=ChemDoodle._AnimatorCanvas.prototype.isRunning;Mol3.dblclick=ChemDoodle.RotatorCanvas.prototype.dblclick;Mol3.timeout=5;Mol3.handle=null; var Fmol='45\ntips\nN -7.117 1.414 1.061\nC -5.883 1.503 1.816\nC -4.825 0.568 1.248\nO -4.335 -0.330 1.933\nC -6.123 1.219 3.298\nC -5.865 2.430 4.184\nC -7.092 3.333 4.198\nC -4.516 3.045 3.826\nO -7.550 3.743 3.108\nO -7.605 3.624 5.295\nO -4.446 4.273 3.622\nO -3.525 2.287 3.734\nH -7.815 0.792 1.350\nH -5.523 2.518 1.717\nH -5.467 0.423 3.612\nH -7.149 0.911 3.435\nH -5.762 2.055 5.191\nH -6.912 3.794 5.937\nH -2.719 2.805 3.797\nN -9.246 3.307 -0.925\nC -8.653 1.995 -0.711\nC -7.319 2.165 -0.007\nO -6.489 2.980 -0.425\nC -8.442 1.253 -2.032\nC -7.601 0.001 -1.854\nO -8.029 -0.906 -1.107\nO -6.502 -0.065 -2.443\nH -8.662 4.116 -0.983\nH -9.317 1.423 -0.080\nH -9.401 0.969 -2.438\nH -7.939 1.908 -2.731\nH -10.238 3.399 -1.014\nH -5.791 0.007 -1.802\nN -4.475 0.798 -0.010\nC -3.472 0.008 -0.700\nC -2.121 0.196 -0.002\nO -1.233 -0.690 -0.101\nC -3.408 0.469 -2.164\nS -2.963 -0.811 -3.398\nH -4.906 1.534 -0.497\nH -3.762 -1.032 -0.658\nH -4.376 0.853 -2.443\nH -2.686 1.266 -2.238\nH -1.932 1.081 0.569\nH -1.659 -0.900 -3.484\n'; Mol3.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol3.startAnimation();Mol3.stopAnimation();function setProj3(yesPers){Mol3.specs.projectionPerspective_3D=yesPers;Mol3.setupScene();Mol3.repaint()}function setModel3(model){Mol3.specs.set3DRepresentation(model);Mol3.setupScene();Mol3.repaint()}function setSpeed3(){Mol3.timeout=500-document.getElementById('spd3').value;Mol3.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol3.startAnimation()}</script> 视图: <input type='radio' name='group2' onclick='setProj3(true)'>投影 <input type='radio' name='group2' onclick='setProj3(false)' checked=''>正交速度: <input type='range' id='spd3' min='1' max='500' onchange='setSpeed3()'/> 模型: <input type='radio' name='model' onclick='setModel3('Ball and Stick')' checked=''>球棍 <input type='radio' name='model' onclick='setModel3('van der Waals Spheres')'>范德华球 <input type='radio' name='model' onclick='setModel3('Stick')'>棍状 <input type='radio' name='model' onclick='setModel3('Wireframe')'>线框 <input type='radio' name='model' onclick='setModel3('Line')'>线型 <input type='checkbox' onclick='Mol3.specs.atoms_displayLabels_3D=this.checked;Mol3.repaint()'>名称 左键: 转动滚轮: 缩放双击: 自动旋转开关 Alt+左键: 移动 <figurecaption>Fig.3</figurecaption></figure>

第四步获取残基的拓扑文件

使用与创建AMBER GAFF力场拓扑文件类似的方法获取残基分子的拓扑文件, 详细过程可参考使用AmberTools+ACPYPE+Gaussian创建小分子GAFF力场的拓扑文件.

作为示例, 我们使用AM1-BCC电荷. 下面的脚本将上面博文中的几个步骤进行了整合, 方便使用.

<table class="highlighttable"><th colspan="2" style="text-align:left">cgu.bsh</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">mol=cgu_h # pdb文件名称 # 1. 判断原子类型, 计算电荷. # 注意我们要使用AMBER而不是GAFF的原子类型 # 我们没有指定体系的电荷和多重度, 因为它们与默认的电荷0, 多重度1一致 antechamber -i $mol.pdb -fi pdb -o $mol.mol2 -fo mol2 -c bcc -at amber -pf y # 2. 检查补充缺失参数 parmchk2 -i $mol.mol2 -f mol2 -o $mol.mod # 3. 生成AMBER参数文件和坐标文件 echo " source leaprc.ff99SBildn loadamberparams $mol.mod mol=loadmol2 $mol.mol2 check mol saveamberparm mol $mol.prm $mol.crd quit" | sleap # 4. 生成GROMACS拓扑文件和坐标文件 acpype -p $mol.prm -x $mol.crd -d </pre></div> </td></tr></table>

运行完成后, 注意检查AMBER优化后的构型sqm.pdb是否正常. 构型正常情况下后续步骤得到的结果才正确. 否则的话, 需要找出原因并修正错误.

完成本步骤后, 我们会得到残基分子的拓扑文件cgu_h_GMX.top.

我们优化后的构型如下, 看起来没有问题.

<figure><script>var Mol4=new ChemDoodle.TransformCanvas3D('Mol-4',642,396);Mol4.specs.shapes_color='#fff';Mol4.specs.backgroundColor='black';Mol4.specs.set3DRepresentation('Ball and Stick');Mol4.specs.projectionPerspective_3D=false;Mol4.specs.compass_display=true; /*//Mol4.specs.atoms_resolution_3D=15; //Mol4.specs.bonds_resolution_3D=15; //Mol4.specs.crystals_unitCellLineWidth=1.5;*/ Mol4.nextFrame=function(delta){var matrix=[];ChemDoodle.lib.mat4.identity(matrix);var change=delta*Math.PI/15000;ChemDoodle.lib.mat4.rotate(matrix,change,[1,0,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,1,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,0,1]);ChemDoodle.lib.mat4.multiply(this.rotationMatrix, matrix)}; Mol4.startAnimation=ChemDoodle._AnimatorCanvas.prototype.startAnimation;Mol4.stopAnimation=ChemDoodle._AnimatorCanvas.prototype.stopAnimation;Mol4.isRunning=ChemDoodle._AnimatorCanvas.prototype.isRunning;Mol4.dblclick=ChemDoodle.RotatorCanvas.prototype.dblclick;Mol4.timeout=5;Mol4.handle=null; var Fmol='45\ntips\nN -8.519 4.036 -0.712\nC -8.148 2.658 -0.500\nC -6.757 2.578 0.170\nO -5.798 3.258 -0.255\nC -8.074 1.787 -1.770\nC -7.923 0.343 -1.421\nO -8.245 -0.247 -0.385\nO -7.343 -0.401 -2.403\nH -7.763 4.565 -1.102\nH -8.938 2.210 0.179\nH -9.021 1.893 -2.365\nH -7.218 2.110 -2.421\nH -9.325 4.099 -1.300\nH -7.296 -1.333 -2.128\nN -6.620 1.715 1.231\nC -5.364 1.547 1.914\nC -4.446 0.508 1.209\nO -4.086 -0.536 1.784\nC -5.523 1.096 3.373\nC -6.229 2.093 4.279\nC -7.726 2.092 4.055\nC -5.617 3.463 4.242\nO -8.346 1.173 3.514\nO -8.527 3.077 4.530\nO -6.156 4.549 4.487\nO -4.287 3.492 3.965\nH -7.367 1.092 1.468\nH -4.813 2.540 1.897\nH -4.487 0.922 3.775\nH -6.075 0.118 3.411\nH -6.086 1.729 5.346\nH -8.019 3.827 4.881\nH -3.969 4.412 3.987\nN -4.040 0.819 -0.069\nC -3.251 -0.089 -0.867\nC -1.777 0.324 -0.887\nO -1.085 0.265 -1.902\nC -3.798 -0.242 -2.291\nS -3.394 -1.804 -3.005\nH -4.308 1.685 -0.480\nH -3.274 -1.111 -0.364\nH -4.920 -0.219 -2.296\nH -3.426 0.595 -2.930\nH -1.363 0.653 0.091\nH -2.104 -1.613 -3.249\n'; Mol4.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol4.startAnimation();Mol4.stopAnimation();function setProj4(yesPers){Mol4.specs.projectionPerspective_3D=yesPers;Mol4.setupScene();Mol4.repaint()}function setModel4(model){Mol4.specs.set3DRepresentation(model);Mol4.setupScene();Mol4.repaint()}function setSpeed4(){Mol4.timeout=500-document.getElementById('spd4').value;Mol4.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol4.startAnimation()}</script> 视图: <input type='radio' name='group2' onclick='setProj4(true)'>投影 <input type='radio' name='group2' onclick='setProj4(false)' checked=''>正交速度: <input type='range' id='spd4' min='1' max='500' onchange='setSpeed4()'/> 模型: <input type='radio' name='model' onclick='setModel4('Ball and Stick')' checked=''>球棍 <input type='radio' name='model' onclick='setModel4('van der Waals Spheres')'>范德华球 <input type='radio' name='model' onclick='setModel4('Stick')'>棍状 <input type='radio' name='model' onclick='setModel4('Wireframe')'>线框 <input type='radio' name='model' onclick='setModel4('Line')'>线型 <input type='checkbox' onclick='Mol4.specs.atoms_displayLabels_3D=this.checked;Mol4.repaint()'>名称 左键: 转动滚轮: 缩放双击: 自动旋转开关 Alt+左键: 移动 <figurecaption>Fig.4</figurecaption></figure>

第五步调整残基的电荷

我们需要明确, 对GROMACS残基数据库中的每个非离子状态的残基, 其净电荷应当为零, 这样它们组合成大分子后净电荷才正确. 但我们在上一步中计算出来的残基的净电荷, 除去用于封端的相邻残基后, 剩余部分的净电荷很可能不是零, 这就需要我们进行处理. 处理方法分两种情况:

如果采用RESP拟合电荷, 可以使用约束, 将每个残基的净电荷约束为零.
如果采用AM1-BCC电荷, 简单的处理方法是将相邻残基的净电荷加到相应的连接原子上, 或者将残基中所有原子的电荷减去其整体净电荷的其平均值.

这里我们使用调整连接原子电荷的方法. 所得拓扑文件中涉及电荷部分如下

<table class="highlighttable"><th colspan="2" style="text-align:left">cgu_h_GMX.top</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">[ moleculetype ] ;name nrexcl cgu_h 3 [ atoms ] ; nr type resi res atom cgnr charge mass ; qtot bond_type 1 N 1 MOL N 1 -0.558901 14.01000 ; qtot -0.559 2 CT 1 MOL CA 2 0.038700 12.01000 ; qtot -0.520 3 C 1 MOL C 3 0.642101 12.01000 ; qtot 0.122 4 O 1 MOL O 4 -0.602101 16.00000 ; qtot -0.480 5 CT 1 MOL CB 5 -0.069400 12.01000 ; qtot -0.550 6 CT 1 MOL CG 6 -0.202700 12.01000 ; qtot -0.752 7 C 1 MOL CD1 7 0.641601 12.01000 ; qtot -0.111 8 C 1 MOL CD2 8 0.641601 12.01000 ; qtot 0.531 9 O 1 MOL OE11 9 -0.524501 16.00000 ; qtot 0.006 10 OH 1 MOL OE12 10 -0.579101 16.00000 ; qtot -0.573 11 O 1 MOL OE21 11 -0.524501 16.00000 ; qtot -1.097 12 OH 1 MOL OE22 12 -0.579101 16.00000 ; qtot -1.676 13 H 1 MOL H 13 0.362500 1.00800 ; qtot -1.314 14 H1 1 MOL HA 14 0.098700 1.00800 ; qtot -1.215 15 HC 1 MOL HB1 15 0.098200 1.00800 ; qtot -1.117 16 HC 1 MOL HB2 16 0.098200 1.00800 ; qtot -1.019 17 HC 1 MOL HG 17 0.121700 1.00800 ; qtot -0.897 18 HO 1 MOL HE1 18 0.460000 1.00800 ; qtot -0.437 19 HO 1 MOL HE2 19 0.460000 1.00800 ; qtot 0.023 ; 以上为CGU 20 NT 1 MOL N_ 20 -0.900801 14.01000 ; qtot -0.878 21 CT 1 MOL C1 21 0.120500 12.01000 ; qtot -0.757 22 C 1 MOL C_ 22 0.660101 12.01000 ; qtot -0.097 23 O 1 MOL O1 23 -0.637101 16.00000 ; qtot -0.734 24 CT 1 MOL C3 24 -0.178400 12.01000 ; qtot -0.913 25 C 1 MOL C5 25 0.633101 12.01000 ; qtot -0.280 26 O 1 MOL OD1 26 -0.546001 16.00000 ; qtot -0.826 27 OH 1 MOL OD2 27 -0.604101 16.00000 ; qtot -1.430 28 H 1 MOL H1 28 0.370300 1.00800 ; qtot -1.059 29 H1 1 MOL H7 29 0.104700 1.00800 ; qtot -0.955 30 HC 1 MOL H9 30 0.094200 1.00800 ; qtot -0.861 31 HC 1 MOL HB3 31 0.094200 1.00800 ; qtot -0.766 32 H 1 MOL H2 32 0.370300 1.00800 ; qtot -0.396 33 HO 1 MOL H3 33 0.455000 1.00800 ; qtot 0.059 ; 以上为3ASP 34 N 1 MOL N__ 34 -0.568901 14.01000 ; qtot -0.510 35 CT 1 MOL C2 35 -0.020300 12.01000 ; qtot -0.530 36 C 1 MOL C__ 36 0.554901 12.01000 ; qtot 0.025 37 O 1 MOL O2 37 -0.530101 16.00000 ; qtot -0.505 38 CT 1 MOL C4 38 -0.010300 12.01000 ; qtot -0.516 39 SH 1 MOL SG 39 -0.369900 32.06000 ; qtot -0.886 40 H 1 MOL H4 40 0.340500 1.00800 ; qtot -0.545 41 H1 1 MOL H8 41 0.124700 1.00800 ; qtot -0.420 42 H1 1 MOL H10 42 0.088200 1.00800 ; qtot -0.332 43 H1 1 MOL H11 43 0.088200 1.00800 ; qtot -0.244 44 HA 1 MOL H5 44 0.016200 1.00800 ; qtot -0.228 45 HS 1 MOL H6 45 0.227800 1.00800 ; qtot 0.000 ; 以上为5CYS </pre></div> </td></tr></table>

检查拓扑文件可以看到, 实际上分子中每个残基的净电荷都很接近零. 其中CGU残基的净电荷为0.023, 很接近零, 不调整的话也影响不大, 但我们还是示例下如何调整.

调整电荷的方法是: 分别计算与CGU残基N端和C端相邻的残基的净电荷, 然后将其值分别累加到CGU残基的N原子和C原子上. 具体来说, 与CGU-N1相邻的是3ASP残基, 将其包含的每个原子的电荷累加, 得到其净电荷为0.0359980, 所以CGU-N1的电荷应当加上这个值, 调整为-0.558901+0.0359980=-0.522903. 同理, CGU-C3的电荷应调整为0.642101-0.0590010=0.5831. 这样调整之后CGU的净电荷就为零了.

第六步整理残基的rtp条目

调整电荷后, 需要将拓扑文件整理为rtp格式, 这样才能将残基的拓扑数据添加到GROMACS力场的残基数据库中. rtp文件的详细说明见GROMACS中文手册第五章 5.6.1 残基数据库.

简单来说, rtp中需要定义的是残基中每个原子的名称, 类型和电荷, 以及原子间的成键相互作用项, 如原子间的连接信息, 键角信息, 二面角信息. 其实, 我们只要定义原子间的连接信息即可, 因为GROMACS可以根据原子间的连接信息自动推断出键角和二面角信息. 但有些特殊的成键相互作用, 如异常二面角, 与相邻残基形成的二面角仍然需要明确定义. 此外, 在定义这些成键相互作用项时可以忽略参数, 但考虑到ambertools给出的参数可能与GROMACS自带力场中的参数存在差异(特别是使用不同版本amber力场时), 再加上有些参数在GROMACS自带的力场中可能并不存在, 所以还是建议将ambertools给出的参数保留在rtp中.

首先将cgu_h_GMX.top另存为cgu.rtp.

删除不需要的内容, 只保留下列部分的内容: [ atoms ], [ bonds ], [ angles ], [ dihedrals ] ; propers, [ dihedrals ] ; impropers.

先整理[ atoms ]部分. 在rtp中, 这一部分只需要保留属于CGU残基的原子即可, 编号为1-19, 因为我们在前面的步骤中已经调整好顺序了, 此外也可以通过N_标识明显看出. rtp中只需要四列: 原子名称, 原子类型, 电荷, 编号. 这些信息在原来的拓扑文件中都是有的, 只不过顺序不同而已, 很容易使用列编辑功能整理好. 注意, CGU残基N1和C3原子的电荷应当是调整过的.

<table class="highlighttable"><th colspan="2" style="text-align:left">cgu.rtp</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">[ atoms ] ; 名称类型电荷编号 N N -0.522903 1 CA CT 0.038700 2 C C 0.5831 3 O O -0.602101 4 CB CT -0.069400 5 CG CT -0.202700 6 CD1 C 0.641601 7 CD2 C 0.641601 8 OE11 O -0.524501 9 OE12 OH -0.579101 10 OE21 O -0.524501 11 OE22 OH -0.579101 12 H H 0.362500 13 HA H1 0.098700 14 HB1 HC 0.098200 15 HB2 HC 0.098200 16 HG HC 0.121700 17 HE1 HO 0.460000 18 HE2 HO 0.460000 19 </pre></div> </td></tr></table>

整理[ bonds ]部分. rtp文件中这一部分使用原子名称指定连接信息, 而不是像原拓扑文件中使用原子编号指定连接信息. 而我们得到的拓扑文件中最后的注释部分就是以原子名称指定的, 所以我们只要将这些内容提取出来, 稍加整理即可.

由于我们只需要保留与CGU残基有关的信息, 所以成键的两个原子中如果有任何一个编号大于19我们就将其删除. 根据原子编号的排列规律, 容易看出第二个的原子编号总是大于第一个, 所以我们直接删除第一个原子编号大于19的左右行即可.

将原拓扑文件中注释部分提取出来放到最前面, 删除原子名称中间的-字符, 同时删除标识势能函数类型的数字1.

最后, 将处于CGU前端并与其相邻的ASP残基的碳原子C_的名称改为-C, 指明是CGU的N原子与前一氨基酸残基的C原子相连. 值得注意的是, 并不需要为CGU残基末端的C原子定义C +N来指定与下一相邻残基的连接. 因为rtp中只定义每个残基的左边(N端)或者开始位置的连接, 末端连接是由下一个残基定义的. 所以删除C N__一行.

<table class="highlighttable"><th colspan="2" style="text-align:left">cgu.rtp</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">[ bonds ] N CA 1.4600e-01 2.7665e+05 N H 1.0090e-01 3.4326e+05 N -C 1.3450e-01 4.0016e+05 CA C 1.5080e-01 2.7472e+05 CA CB 1.5350e-01 2.5363e+05 CA HA 1.0930e-01 2.8108e+05 C O 1.2140e-01 5.4225e+05 CB CG 1.5350e-01 2.5363e+05 CB HB1 1.0920e-01 2.8225e+05 CB HB2 1.0920e-01 2.8225e+05 CG CD1 1.5080e-01 2.7472e+05 CG CD2 1.5080e-01 2.7472e+05 CG HG 1.0920e-01 2.8225e+05 CD1 OE11 1.2140e-01 5.4225e+05 CD1 OE12 1.3060e-01 3.9028e+05 CD2 OE21 1.2140e-01 5.4225e+05 CD2 OE22 1.3060e-01 3.9028e+05 OE12 HE1 9.7400e-02 3.0928e+05 OE22 HE2 9.7400e-02 3.0928e+05 </pre></div> </td></tr></table>

对于[ angles ]和[ dihedrals ] ; propers部分, 处理方法类似[ bonds ]部分. 直接删除包含非CGU原子的项即可, 而且无需考虑与前后残基的连接原子. 有时我们甚至可以直接将这两部分删除, 但这样做有可能出现相互作用函数缺失的问题, 所以建议还是尽量保留.

<table class="highlighttable"><th colspan="2" style="text-align:left">cgu.rtp</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">[ angles ] N CA C 1.1156e+02 5.5789e+02 N CA CB 1.1213e+02 5.5103e+02 N CA HA 1.0932e+02 4.1689e+02 CA N H 1.1678e+02 3.8526e+02 CA C O 1.2311e+02 5.6928e+02 CA CB CG 1.1063e+02 5.2894e+02 CA CB HB1 1.1005e+02 3.8802e+02 CA CB HB2 1.1005e+02 3.8802e+02 C CA CB 1.1053e+02 5.3379e+02 C CA HA 1.0766e+02 3.9857e+02 CB CA HA 1.1007e+02 3.8794e+02 CB CG CD1 1.1053e+02 5.3379e+02 CB CG CD2 1.1053e+02 5.3379e+02 CB CG HG 1.1005e+02 3.8802e+02 CG CB HB1 1.1005e+02 3.8802e+02 CG CB HB2 1.1005e+02 3.8802e+02 CG CD1 OE11 1.2311e+02 5.6928e+02 CG CD1 OE12 1.1220e+02 5.8442e+02 CG CD2 OE21 1.2311e+02 5.6928e+02 CG CD2 OE22 1.1220e+02 5.8442e+02 CD1 CG CD2 1.1161e+02 5.3605e+02 CD1 CG HG 1.0968e+02 3.9497e+02 CD1 OE12 HE1 1.0737e+02 4.2836e+02 CD2 CG HG 1.0968e+02 3.9497e+02 CD2 OE22 HE2 1.0737e+02 4.2836e+02 OE11 CD1 OE12 1.2288e+02 6.4752e+02 OE21 CD2 OE22 1.2288e+02 6.4752e+02 HB1 CB HB2 1.0835e+02 3.2995e+02 [ dihedrals ] ; propers N CA C O 180.00 0.00000 2 N CA CB CG 0.00 0.65084 3 N CA CB HB1 0.00 0.65084 3 N CA CB HB2 0.00 0.65084 3 CA CB CG CD1 0.00 0.65084 3 CA CB CG CD2 0.00 0.65084 3 CA CB CG HG 0.00 0.66944 3 C CA CB CG 0.00 0.65084 3 C CA CB HB1 0.00 0.65084 3 C CA CB HB2 0.00 0.65084 3 CB CA C O 180.00 0.00000 2 CB CG CD1 OE11 180.00 0.00000 2 CB CG CD1 OE12 180.00 0.00000 2 CB CG CD2 OE21 180.00 0.00000 2 CB CG CD2 OE22 180.00 0.00000 2 CG CD1 OE12 HE1 180.00 9.62320 2 CG CD2 OE22 HE2 180.00 9.62320 2 CD1 CG CD2 OE21 180.00 0.00000 2 CD1 CG CD2 OE22 180.00 0.00000 2 CD2 CG CD1 OE11 180.00 0.00000 2 CD2 CG CD1 OE12 180.00 0.00000 2 OE11 CD1 OE12 HE1 0.00 7.94960 1 OE11 CD1 OE12 HE1 180.00 9.62320 2 OE21 CD2 OE22 HE2 0.00 7.94960 1 OE21 CD2 OE22 HE2 180.00 9.62320 2 H N CA C 0.00 0.00000 0 H N CA CB 0.00 0.00000 0 H N CA HA 0.00 0.00000 0 HA CA C O 0.00 0.00000 0 HA CA C O 0.00 3.34720 1 HA CA C O 180.00 0.33472 3 HA CA CB CG 0.00 0.65084 3 HA CA CB HB1 0.00 0.65084 3 HA CA CB HB2 0.00 0.65084 3 HB1 CB CG CD1 0.00 0.65084 3 HB1 CB CG CD2 0.00 0.65084 3 HB1 CB CG HG 0.00 0.62760 3 HB2 CB CG CD1 0.00 0.65084 3 HB2 CB CG CD2 0.00 0.65084 3 HB2 CB CG HG 0.00 0.62760 3 HG CG CD1 OE11 0.00 0.00000 0 HG CG CD1 OE11 0.00 3.34720 1 HG CG CD1 OE11 180.00 0.33472 3 HG CG CD1 OE12 180.00 0.00000 2 HG CG CD2 OE21 0.00 0.00000 0 HG CG CD2 OE21 0.00 3.34720 1 HG CG CD2 OE21 180.00 0.33472 3 HG CG CD2 OE22 180.00 0.00000 2 </pre></div> </td></tr></table>

整理[ impropers ]部分. 首先将[ dihedrals ] ; impropers修改为[impropers] , 然后删除存在编号大于19的项, 但保留C_- CA- N- H和CA- N__- C- O对应的项. 这样共保留4项, 其中两项对应于两个羧基上的氢原子, 剩下两项对应于CGU与前后残基所成肽键的异常二面角项. 同样, 我们需要将后面两项中的C_改为-C, N__改为+N指明是前后残基的相连原子, 而不是CGU自身的原子.

<table class="highlighttable"><th colspan="2" style="text-align:left">ahk</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%">1 2 3 4 5</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">[ impropers ] CA +N C O 180.00 4.60240 2 CG OE11 CD1 OE12 180.00 4.60240 2 CG OE21 CD2 OE22 180.00 4.60240 2 -C CA N H 180.00 4.60240 2 </pre></div> </td></tr></table>

成键部分整理完成后, 在文件的开头添加默认成键类型, CGU残基的定义, 最终得到了满足要求的rtp文件, 其整个结构如下

<table class="highlighttable"><th colspan="2" style="text-align:left">cgu.rtp</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">[ bondedtypes ] 1 1 9 4 1 3 1 0 [ CGU ] [ atoms ] ... [ bonds ] ... [ angles ] ... [ dihedrals ] ; propers ... [ impropers ] ... </pre></div> </td></tr></table>

第七步准备残基的hdb文件

我们在使用pdb2gmx的时候, 经常需要使用-ignh选项来自动添加氢原子, 所以, 我们也需要为新残基CGU准备一个hdb文件, 用于说明如何自动添加氢原子. GROMACS手册第五章 5.6.4 氢数据库对此有说明, 但不够直观. 下面给出氢原子添加方式的图形化说明.

创建cgu.hdb文件, 根据CGU中氢原子的名称

定义hdb文件如下:

其中第一行为残基名称, 待添加氢的类型数, 接下来每行列出一类待添加氢的信息: 加氢数目, 加氢方式, 加氢名称, 原子i, 原子j, 原子k, 原子l.

第八步模拟测试

我们使用的是amber99sb-ildn力场, 所以将准备好的cgu.rtp和cgu.hdb文件放到 /GMX安装路径/share/gromacs/top/amber99sb-ildn.ff目录下, 这样CGU残基的信息就被添加到力场数据库中了. 这里要说明的是, 我们不建议直接将rtp和hdb里面的信息添加到GROMACS力场自带的aminoacids.rtp和aminoacids.hdb文件中, 那样做容易污染原先的文件.

需要注意的是, 由于原始gid.pdb中也含有HYP残基, 虽然我们使用的amber99sb-ildn力场中含有它的拓扑数据, 但因为它不是常见的标准残基, 所以我们要查看一下它的原子名称是否与残基数据库中的一致. 打开/GMX安装路径/share/gromacs/top/amber99sb-ildn.ff/aminoacids.rtp, 搜索HYP就可以看到HYP残基中每个原子名称, 与gid.pdb中的名称进行对比, 可以发现, gid.pdb中HYP的CD必须改为CD2才与力场数据库一致.

新残基添加到数据库中以后, 我们就可以使用pdb2gmx来测试是否添加成功了. 运行

选择amber99sb-ildn力场, TIP3P水模型. 如果不能通过, 仔细检查cgu.rtp, 修正错误, 直至成功通过, 得到conf.gro和topol.top文件.

使用可视化软件打开conf.gro文件, 检查自动添加的氢原子是否正确. 如果不正确, 需要修改cgu.hdb, 重复上面的pdb2gmx命令, 直至添加的氢原子正确为止.

按照蛋白模拟的常规流程进行短时间的模拟, 检查轨迹是否有明显的异常(如原子跳跃异常, 成键断开等). 如果有明显异常, 寻找产生问题的原因. 修正后, 重复模拟, 直至轨迹正常.

其他可能的问题

我们在处理过程中, 使用了amber的原子类型, 目的在于尽量避免引入新的原子类型. 如果确实需要引入新的原子类型, 我们需要添加新的原子类型及其参数到力场文件中, 还需要定义新原子类型与其他原子类型的成键相互作用参数. 如果混用不同的力场, 很多成键参数是没有的, 处理起来很麻烦. 因此建议尽量避免引入新的原子类型, 不要混用力场.
如果pdb2gmx执行成功, 而grompp时出错, 提示存在未定义的成键相互作用, 那么必须定义这些缺失的成键相互作用参数. 这些参数或者来自文献, 或者自行拟合, 或者使用力场中的类似参数. 通常少数几个成键参数的不同对结构和分子间相互作用的影响不大. 但情况并非总是这样, 所以你还是要留意这一点.

仓库源文，站点原文

第一步 获取蛋白的PDB文件

第二步 更新残基类型文件residuetypes.dat

第三步 创建包含残基的完整分子

第四步 获取残基的拓扑文件

第五步 调整残基的电荷

第六步 整理残基的rtp条目

第七步 准备残基的hdb文件

第八步 模拟测试