本教程用于引导学生运行蛋白质的分子动力学(MD)模拟.如果待研究蛋白不包含非标准氨基酸残基,教程所用的流程可作为研究蛋白的一个良好起点.完成本教程后,学生应该能够知晓设置和运行模拟所涉及的步骤,以及在不同模拟阶段应如何选择.此外,学生还应该知道如何对模拟结果进行核查以保证模拟质量,并了解可以使用哪些分析方法,哪些分析方法更适合自己的需要.
本教程的目的在于考察一段小肽结合到其蛋白受体时,生物分子相互识别的两种机理–构象选择以及诱导匹配–所具有的不同贡献.完成本教程后,学生应该能够:
Linux下的所有命令都需要小心输入,因为shell(负责解析命令的程序)是区分大小写的.另一经常发生的错误是将数字0当成了大写字母O,小写字母l当成了数字1,或者反过来.你可以简单地复制粘贴这些命令:先使用鼠标选中它们,然后在需要输入的地方按下鼠标中键(Linux模式,Windows下未必可行).首先来试试下面的命令:
ls-l上面的命令列出当前工作目录下的所有文件和文件夹.当遇到filenotfound之类的错误时,可以使用这个命令来检查.
pwd这个命令给出当前目录的完整路径.
cdDesktop此命令展示了如何更换目录.
在决定细胞命运的大量蛋白-蛋白相互作用中,多肽起着关键作用,并占了大约40%的比例.从辅激活剂到抑制剂,它们在许多信号和识别途径中都有涉及,并且已经被证实与大量蛋白结构域存在相互作用,例如,MHC,SH3和PDZ结构域就是由于其多肽亲和能力而知名.由于功能的多样性以及间接参与的许多生物途径的重要性,使得多它们与许多疾病相联系.药物设计中的一个领域就专门研究多肽并用于治疗各种疾病.与小分子抑制剂相比,多肽的优点在于可以模拟蛋白结合的结构域,并且自身足够大能竞争性地抑制蛋白蛋白相互作用.许多药物先导分子就包含抗菌多肽.
尽管科学家们收集了蛋白多肽相互作用的大量数据,确定它们形成的复合物的结构仍然是个挑战.这主要源于两个障碍,多肽柔性极高,并且与其底物的相互作用很弱,这突出了它们在信号传导或调节中的重要性,因为这些功能通常依赖于瞬态过程.这些障碍使得实验结构测定并不简单,需要借助额外的计算方法如生物分子对接才能完成.从建模的角度看,常规的用于蛋白-配体或者蛋白-蛋白对接的算法也经常受到柔性问题的困扰.
解决柔性问题的一个途径是进行分子动力学模拟,对多肽的构象空间进行采样,收集具有代表性的构象,并用于预测它们与蛋白受体的相互作用.
上个世纪人们提出了一些理论来解释分子的识别过程.在这些理论中,构象选择和诱导匹配机理在过去的50年中得到了大致相同多的支持.构象选择假定,在配体不存在的情况下,蛋白处于多种离散构象状态的平衡中,其中包括倾向与配体结合的那个构象.这一概念与诱导匹配理论相反,后一理论起初引入用于描述酶的活性,认为构象匹配是由与底物的结合导致的.
在本教程中,我们要考察这一两难问题–诱导匹配和/或构象选择,方法是使用不同初始构象的分子动力学模拟对多肽的构象空间进行采样.我们将不同的轨迹与实验确定的多肽在与其蛋白受体形成的复合物中的构象(结合构象)进行比较,并检查这一构象在MD过程中是否出现,这样就能验证构象选择理论.
默认情况下,你的系统已经正确设置好了.
你可以使用gmxluck命令测试GROMACS是否正确安装.幸运的话,你会得到一句引言.这说明你的GROMACS已经正确安装好了.
为开始教程,根据多肽的PDB号下载相应的结构,并记下下列信息:
首先使用分子查看软件来看看分子的结构,我们建议使用PyMOL.可使用下面的PyMOL命令来载入结构
PyMOL1klu.pdb执行上面的命令后,PyMOL会启动,在一个窗口中使用线形模式来显示分子,并在主窗口的右面列出各个对象,这些对象可以通过点击其名称来移除.每个对象名称的后面是菜单,用于更改显示模式.试着使用卡通模式来显示结构.对那些习惯使用键盘的人(我们强烈建议你也使用这种方式),可以通过在窗口中键入下面的命令来达到目的
showcartoon使用上一步骤中收集到的信息,通过键入下面的命令来选择包含多肽的pdb链
selectsel_peptide,chainX然后抽取多肽链保存到另一个对象
extractpeptide,sel_peptide计算并显示蛋白的溶剂可及表面
assurface,1klu对蛋白表面和多肽使用两种可区分的颜色进行着色
colorwhite,1klucolortv_blue,peptide选择处于多肽5埃之内的所有蛋白残基进行着色,这构成了蛋白-多肽的界面
selectinterface,byres(1kluwithin5Aofpeptide)colorred,interface如果对显示结果满意,你可以进一步使用ray命令改善显示结果,并使用pngfilename命令保存图片.
bg_colorwhitesetopaque_background,offray1000pngmy_pretty_peptide.png使用下面的命令记下在PyMOL中显示的结合多肽的氨基酸序列
setseq_view,on比较PyMOL显示的序列和蛋白数据库给出的序列,二者是否相同为什么
PyMOL也可以用于生成多肽的初始构象.默认情况下,你可以使用Build菜单,选择构象和需要增加的残基.此外,我们提供了一个脚本build_seq,可以对给定的氨基酸序列生成相应的理想构象.这个脚本可能很有用,因为对于一些多肽构象,如聚脯氨酸II,PyMOL并不知道所需要的phi/psi角度组合.你可以通过在窗口中键入下面的命令来获知如何使用这个脚本
helpbuild_seq使用这个脚本生成多肽的三个不同构象,我们在教程中要使用它们,所以要给它们指定不同的名字,将下面命令中的XXXXXX替换为结合多肽的序列(PyMOL中有显示)
build_seqpeptide_helix,XXXXXX,helixbuild_seqpeptide_extended,XXXXXX,betabuild_seqpeptide_polyproline,XXXXXX,polypro使用卡通模式显示新创建的多肽
showcartoon,peptide_*使用align命令在空间中对齐并叠合不同的构象,查看它们的差异程度有多大
alignpeptide_helix,peptidealignpeptide_extended,peptidealignpeptide_polyproline,peptide对齐与叠合计算会返回原子位置的根均方偏差(RMSD),它定量地表征了不同构象的结构异同.请记下每个align命令给出的RMSD值.
最后,使用save命令保存构象,类似下面保存螺旋多肽的命令
savehelical_peptide.pdb,peptide_helix得到的pdb文件可以用于开始准备模拟.
现在使用命令quit退出PyMOL.正如你注意到的,绘制结构必需的所有信息都保存在.pdb文件中.你可以使用gedit来查看pdb文件,并试着理解这种文件格式.
pdb文件包含了很多涉及蛋白,用于确定结构的实验方法和条件等的信息.PDB数据库网页实际就是将这些信息展示给你的.pdb文件中也包含了每一个原子的直角坐标.注意,pdb文件中一般不包含成键连接信息,而PyMOL像大多数分子查看软件一样绘制了原子之间的成键,这些成键是自动根据原子间的距离来确定的.
结合态多肽的pdb文件和你自己使用build_seq生成的有什么不同(如氢原子)
为你的结构创建一个目录.由于最后需要将结果合并起来,最好通过组合PDB标识(如peptide_bound,peptide_helix等)和个人标识指定明确且唯一的目录名称.将结构文件复制到目录中并切换目录.举例来说,
mkdirpeptide_boundcpbound_peptide.pdbpeptide_bound/cdpeptide_bound/现在开始真正的MD部分.记住在每一步都要指定正确的文件名.特别注意,教程简单地使用protein.pdb和protein-EM-solvated.gro作为通用名称,这些通用名称应该替换为待研究蛋白的名称(如bound_peptide.pdb和bound_peptide-EM-solvated.gro).
请仔细阅读教程,并检查在每一步是否成功完成.请细心阅读程序输出!!!如果程序给出错误信息,这些信息通常是自明,容易理解.检查文件类型和程序输出以理解每一步的过程.大多数文件是可读的,除了扩展名为.tpr,.xtc,.trr和.edr的文件
准备模拟体系是MD中最重要的步骤.执行MD模拟可以在原子尺度洞察动力学的过程,用于理解实验观察到的现象,验证理论假说,或者为一个有待实验验证的新假说提供基础.然而,对于上述各种情形,设计的模拟必须适合目的,要根据实际情况对模拟过程进行设计,这意味着设置模拟的时候必须十分小心.
缺失残基/侧链/原子,非标准基团
本教程我们模拟的是多肽.整个过程的首要步骤是选择初始结构,如前面所说.然后就要检查这个结构是否缺失残基和原子,这些残基和原子模拟时必须考虑,因此必须想办法补充完整.从pdb库下载的文件中会列出缺失的残基和原子.由于本教程不涉及蛋白建模,因此所用的初始结构是完整的,不存在缺失,因此无需进一步处理.
另一个需要注意的问题是,pdb结构文件中可能包含非标准残基,修饰过的残基,溶剂分子或者配体.这些基团没有力场参数可用.如果存在这样的基团,要么除去它们,要么补充它们的力场参数,这牵涉到MD的高级课题.本教程假定结构不含这类基团,只包含天然氨基酸.
最后,一些晶体散射数据的分辨率足够高,能够在密度网格中区分水分子.通常这些水分子只以氧原子表示,在开始准备步骤前必须去除它们,除非出于特殊目的要保留它们.对本教程的多肽不需要考虑这种情况.
幸运的是,大多数的这些”问题”分子在PDB文件中都被标识为杂原子(HETATM),这样可以使用sed命令轻松地移除它们
sed-i-e"/^HETATM/d"protein.pdb结构质量
一个分子可由各个原子的坐标,成键和非键相互作用的描述来定义.由于从PDB文件获得的结构只含有原子坐标,我们首先必须构建拓扑,它给出了原子类型,电荷,成键情况等信息.拓扑对应着某一特定的力场,选择使用何种力场是一个值得仔细考虑的问题.这里我们使用AMBER99sb-ILDN全原子力场,它被广泛用于采样和折叠模拟的研究.
拓扑与结构的匹配很重要,这意味着结构数据也需要根据所用力场进行转换.可以使用pdb2gmx程序转换结构并建立拓扑数据.该程序可以构建分子的拓扑,前提是这些分子由特定的构建单元(如氨基酸)组成.pdb2gmx使用构建单元的数据库进行转换,对于数据库中没有的分子或残基,程序无法识别,转换过程会失败.输入下面的命令进行结构转换,提示时选择AMBER99SB-ILDN力场和TIP3P(TransferableIntermolecularPotential3P)水分子模型.注意-ignh选项,它表明程序会首先移除初始文件中的氢原子,然后再根据相应力场的说明进行重建.使用-ter选项后GROMACS会提示分子的末端变化,合适的选择应该能够正确地表示体系.在本教程中,这取决于多肽是一个更大蛋白的片段(非带电末端)还是一个完整的分子(末端带电).
gmxpdb2gmx-fprotein.pdb-oprotein.gro-pprotein.top-ignh-ter注意,你可以使用下面的命令将gro文件和pdb文件进行互换
gmxeditconf-fXXXXX.gro-oXXXXX.pdb仔细阅读屏幕输出的提示,并检查组氨酸质子化状态所做的选择,注意蛋白质的总电荷数.也要浏览输入结构的文件(protein.pdb)和输出的结构文件(protein.gro,GROMACS格式).注意两种格式的区别.最后浏览拓扑文件及其结构.
记下转换前后原子数目的区别,解释为什么.选择拓扑文件中的一个氨基酸残基,列出其中的每个原子及其类型,电荷.
gmxeditconf-fprotein.gro-oprotein-PBC.gro-btdodecahedron-d1.2看看editconf的输出,注意体积的变化.另外,也看看protein-PBC.gro文件的最后一行(提示:可使用命令tail-n1protein-PBC.gro).在GROMACS格式(.gro)中,最后一行指定单位晶胞的形状,并总是使用三斜矩阵的表示方法,其中前面的三个数字指定对角元素(xx,yy,zz),后面的六个数字指定非对角元素(xy,xz,yx,yz,zx,zy).
新盒子的体积多大
gmxgrompp-v-f01_em_vac_PME.mdp-cprotein-PBC.gro-pprotein.top-oprotein-EM-vacuum.tprgrompp程序将会提示此体系的净电荷不为零,还会输出一些与体系和控制参数有关的其他信息.该程序也会产生一个额外的输出文件(mdout.mdp),里面包含所有控制参数的设置.下一步,运行mdrun
gmxmdrun-v-deffnmprotein-EM-vacuum由于只有3150个原子,能量最小化很快就完成了.使用选项-v会将每步的势能和最大受力打印出来,以便于跟踪能量最小化过程.-deffnm选项设定了所有输出文件的默认文件名,避免了对每个输出文件进行命名,从而减少了需要设置的选项.现在,结构弛豫好了,我们该添加溶剂了.
使用了什么方法进行能量最小化参数文件中指定了多少步,实际使用了多少步进行能量最小化什么原因致使能量最小化在指定步之前停止体系最终的势能为多少使用PyMOL载入能量最小化之前和之后的结构,并进行比较.
为了更好地查看结构的区别,可以将优化后的结构叠合到初始结构上(已经重命名为protein.pdb).使用align命令可以将结构protein-EM-vacuum.pdb与protein.pdb对齐,
gmxsolvate-cpprotein-EM-vacuum.gro-csspc216.gro-pprotein.top-oprotein-water.gro看看protein.top文件的结束部分,其中有添加溶剂的内容,检查添加的溶剂分子的数目.
向体系中添加了多少水分子,对应的体积多大体系中现在有多少原子
用以前学过的命令,将.gro文件转换为.pdb文件.用PyMOL程序载入溶剂化的结构(protein-water.pdb).在PyMOL中,可用下面的命令显示盒子形状
showcell执行命令后,PyMOL会绘制一个线框来显示三斜晶胞盒子的边界.可能看起来不太明显,但这个三斜形状对应于菱形十二面体.另一个值得注意的事情是,并非所有的溶剂都处于盒子内,但都以长方体块状排布.有时,蛋白质也会部分地伸出水外面.
蛋白伸出水面外部为什么不是问题
gmxgrompp-v-f02_em_sol_PME.mdp-cprotein-water.gro-pprotein.top-oprotein-water.tpr接下来,文件protein-water.tpr可作为genion程序的输入文件.我们设定了NA+/CL-(-pnameNA+-nnameCL-)的浓度为0.15M(-conc0.15),并指定必须加入超量的某一特定种类的离子以使体系呈电中性(-neutral).genion程序会询问用离子取代何种分子,选择SOL组.
gmxgenion-sprotein-water.tpr-oprotein-solvated.gro-conc0.15-neutral-pnameNA+-nnameCL-注意:为中和体系电荷而添加的离子的数量.通过将一些水分子替换为离子,体系的拓扑文件protein.top不再正确.你可以修改拓扑文件,减少溶剂分子的数目,同时在文件的最后一部分(molecules)中SOL的后面增加一行,指定添加的NA离子的数目,并再增加另一行指明CL离子的数目.genion也提供了一个-p选项能够自动更新拓扑文件,更方便.
向体系中添加了多少钠离子和氯离子
到此为止,整个模拟系统已经准备好了.在开始成品模拟前,唯一要做的是对体系再次进行受控的弛豫.由于添加溶剂以及替换离子,可能会产生不利的相互作用,例如,原子间的重叠,同种电荷太接近等,因此需要对体系再次进行能量最小化,其步骤与前面相同:先运行grompp,再运行mdrun:
gmxgrompp-v-f02_em_sol_PME.mdp-cprotein-solvated.gro-pprotein.top-oprotein-EM-solvated.tprgmxmdrun-v-deffnmprotein-EM-solvated参数文件中指定了多少步,优化过程进行了多少步体系最终的势能多大
查看上面提供的mdp文件,你使用哪个温度模拟体系选择这一特定温度的原因何在
编辑完参数文件后,使用grompp和mdrun程序来准备和运行模拟:
gmxgrompp-v-f03_nvt_pr1000_PME.mdp-cprotein-EM-solvated.gro-pprotein.top-oprotein-NVT-PR1000.tprgmxmdrun-v-deffnmprotein-NVT-PR1000为什么我们建议你修改随机数种子的值而不是使用自动生成的随机种子(默认值-1)模拟长度是多少ps在拓扑文件中位置限制是如何包含进去的有两个分别耦合到热浴的组,它们是哪些组,每个组中包含哪些粒子
看看总能量,势能和动能的变化会很有趣.在前一步,粒子没有速度,所以没有动能,因此也没有温度.现在,模拟一开始,原子就被赋予速度因此获得了动能.
gmxenergy-fprotein-NVT-PR1000.edr-othermodynamics-NVT-PR1000.xvg运行上面的命令后,会提示你一系列能量项以供选择,所选能量项的结果将会被写入指定的输出文件中.输入与温度,势能,动能和总能量相对应的数字,最后输入0(零),回车.(提示:你可以通过使用管道(|)自动选择某些项,例如echo101112140|gmxenergy...,或echoTemperature0|gmxernergy...)
xmgr-nxythermodynamics-NVT-PR1000.xvg注意:绿色曲线是总能量,黑色是势能,蓝色是温度,红色是动能.你也可以点击曲线并填写Legend窗口以改变图例.别忘了点击accept以接受改变.
温度如何变化势能/动能/总能量如何变化,如何解释
gmxgrompp-v-f04_npt_pr_PME.mdp-cprotein-NVT-PR1000.gro-pprotein.top-oprotein-NPT-PR1000.tprgmxmdrun-v-deffnmprotein-NPT-PR1000查看控制参数文件,找到控制压力耦合的参数.
运行结束后,再次查看一下能量和温度.使用前面的方法抽取它们.数据提取方法同前:
gmxenergy-fprotein-NPT-PR1000.edr-othermodynamics-NPT-PR1000.xvg观察能量图.
现在开始重复相同的平衡模拟,并逐步放开限制
sed-e's/100010001000/100100100/g'posre.itp>tmp.itpmvtmp.itpposre.itpgmxgrompp-v-f04_npt_pr_PME.mdp-cprotein-NPT-PR1000.gro-pprotein.top-oprotein-NPT-PR100.tprgmxmdrun-v-deffnmprotein-NPT-PR100再减小一次
gmxgrompp-v-f05_npt_NOpr_PME.mdp-cprotein-NPT-PR10.gro-pprotein.top-oprotein-NPT-noPR.tprgmxmdrun-v-deffnmprotein-NPT-noPR运行结束后,再次查看一下能量和温度,同时注意观察压力变化.使用前面的方法抽取它们.数据提取方法同前:
gmxenergy-fprotein-NPT-noPR.edr-othermodynamics-NPT-noPR.xvg看看能量,温度和压力的变化曲线.
温度如何变化压力如何变化
最后一个问题与如何从有限数量的粒子系统中提取热力学性质直接有关.粗略地说,热力学指的是大量粒子(例如,几十亿个而不是几千个)的行为.对大量粒子的性质进行平均能够减少波动,相反,仅对少量粒子进行平均,较大的波动无可避免.
由于我们引入了压力耦合,体系的密度会发生变化.从能量文件中提取密度数据,方法如下:
终于到了最后一步.我们已经得到了或多或少平衡好的溶剂化体系,其中包含我们感兴趣的蛋白质,所以该进行成品模拟了.记住,成品模拟并不表示整个模拟都可以用来分析感兴趣的性质.虽然已经消除了初始构象的一些影响,体系还不太可能已经达到平衡状态.在分析阶段,我们会检查模拟的哪一部分可以认为处于平衡状态,适合用于进一步处理分析.但首先需要设定运行参数.这里只需要运行另一个模拟步骤,类似于准备体系的最后一步.然而,这一步也是另一处需要考虑模拟目的的地方,所以应当选择能满足待分析性质的有关控制参数.可考虑以下问题:
为了完成50纳秒的模拟,需要多少步
gmxgrompp-v-f06_md_PME.mdp-cprotein-NPT-noPR.gro-pprotein.top-oprotein_md.tpr虽然运行输入文件是二进制格式,我们还是可以查看其内容.在某些情况下,模拟会发生一些意外,这可能与内部控制和力场参数有关.在这种情况下,查看运行输入文件尤为有用.gmxdump程序可以将运行输入文件转换为可读格式.转换后的内容可能会有很多页,因此建议将其重定向到一个文件中,或使用more或less命令分屏显示文件.输入以下命令查看运行输入文件的内容(注意,在实际操作时,需要把文件名protein_md.tpr替换为运行grompp后所得到的文件名):
gmxdump-sprotein_md.tpr|&less第三部分分子动力学模拟数据的分析(一)模拟结束后就可以进行数据分析了.这是一个重要的过程,包括三个阶段.首先,有必要进行一些标准的检查,以便对模拟质量进行评估.如果评估结果表明模拟良好,就可以对每个模拟进行分析,回答预设的研究问题了.最后,来自不同模拟的结果可以综合起来.
注:文件名应当能反映文件的内容,根据模拟体系的不同而不同.这里我们假定使用默认的文件名,这意味着会有以下文件:
在进行其他分析之前,必须确认模拟正常结束.有许多原因可能导致模拟中断,特别是涉及力场或体系平衡不充分的问题,要检查模拟是否正常结束,可以使用gmxcheck程序
gmxcheck-ftraj.xtc确认模拟运行了50纳秒.
现在是有趣的部分,尽管多数分析都归结为从轨迹文件中提取图像,MD当然最重要是的体系的运动.我们要看看轨迹.
首先使用GROMACS提供的查看器ngmx来看看.虽然该程序的完善程度和视觉效果不及其他查看软件,但它能够根据拓扑文件的信息绘制成键.其他查看软件可能隐含长程键,导致这些键被认为太长而不画出,或者会在非常接近的原子之间画出键.这是对模拟结果分析的一个常见错误.使用ngmx载入拓扑和轨迹文件:
gmxngmx-stopol.tpr-ftraj.xtc看看程序菜单,试试不同的选项.播放动画.观看过程.通过右边的选项控制.右击或左击选择选项来改变查看.
为了可视化轨迹,我们将从轨迹中抽取1000帧(-dt50),并去除水分子(当提示选择时,选择蛋白Protein).此外,我们还要消除跨过盒子边界的跳跃形成连续轨迹(-pbcnojump).要处理这些事情,我们使用瑞士军刀般的GROMACS工具trjconv,它有1001个可能的组合选项.我们使用它输出一个多模型的pdb文件,用于PyMOL可视化.
当提示选择trjconv的输出组时,选择蛋白Protein(组1)从而忽略溶剂与离子.(提示:如前面所述,可以使用echo来传递选项:echo1|trjconv...)
gmxtrjconv-stopol.tpr-ftraj.xtc-oprotein.pdb-pbcnojump-dt50在PyMOL中载入轨迹
pymolprotein.pdb载入所有帧后,播放动画
mplay当播放动画时,所有其他控制仍可以工作.你可以使用鼠标旋转或缩放体系,也可以改变分子的表示模式.
spectrumshowcell如果多肽扩散超过盒子边界,会发生什么
如果一切正常的话,现在你可以看到多肽的扩散,翻转和扭动.但我们对内部运动更感兴趣,而不是整体行为.在PyMOL中,你可以使用命令intra_fit将轨迹中的所有其他帧与第一帧对齐,随后,你可以使用orient命令设定聚焦点在多肽上
intra_fitproteinorient现在所有帧都已经叠合好了,你可以看到多肽的一些部分比另一些部分运动得更剧烈,这种运动的差异在后面会进行定量化分析.
当然,使用卡通模式表示多肽显示效果会更好:
showcartoon上面的命令会将碳骨架表示为粗管状,而无法显示正确的二级结构元素,因为.pdb文件中没有二级结构信息.PyMOL能够计算蛋白的二级结构,但只计算一帧,然后将结果应用到所有的帧.例如,下述命令可以计算第一帧的二级结构:
dss通过指定状态编号,可以改变要计算的帧
dssstate=1000最后,让我们同时查看所有帧,检查多肽的柔性和刚性区域
如果有足够的机时,你可以考虑制作一段不错的动画.你可能已经注意到了,轨迹的噪声很大.这基本上是热噪声,因此只是蛋白正常行为的一部分,但这些噪声对制作好的动画会有影响.我们可以滤除这些高频的运动,只保留更慢更平滑的整体运动.为此,可以使用filter命令:
gmxfilter-stopol.tpr-ftraj.xtc-olfiltered.pdb-fit-nf5现在在PyMOL中载入滤过后的轨迹,设置二级结构(dss),显示二级结构(showcartoon),隐藏碳骨架上的侧链(hidelines,not(namec,n,o)),以你喜欢的颜色进行着色,然后使用orient命令设置好视角.现在可以开始制作动画了:
mpeg_encodemovie.param3.质量保证对轨迹进行了最初的可视化检查后,该对模拟质量进行一些更彻底的检查了.质量保证(QA)包括对一些热力学参数收敛程度的测试,如温度,压力,势能和动能等.更常用的含义,QA试图评价模拟是否达到平衡.通过起始结构和平均结构的均方根偏差(RSMD)也可以对结构的收敛性进行检查.接下来必须检查相邻周期性映像之间没有相互作用,因为这种相互作用会导致一些非物理效应.最后一步QA测试是计算原子的均方根波动,并与晶体学数据的b因子进行比较.
我们首先从能量文件中提取一些热力学数据,包括温度,压力,势能,动能,晶胞体积,密度以及盒子大小等.这些量中的大多数已经在前面的准备步骤阶段检查过了.能量分析使用energy命令进行,该程序读入能量文件,也就是模拟过程中生成的扩展名为.edr的文件.energy命令会提示需要从能量文件中抽取哪些项,并为其生成一个图形.使用下面的命令
xmgrace-nxytemperature.xvg体系的平均温度和热容多少
通过名称引用能量项可实现能量文件的自动处理.使用echo和管道(|)可以将一个程序输出重定向为另一个程序的输入,这样energy的选择可以自动完成.要抽取多个项,各项之间必须以\n分割,复制粘贴或键入以下命令行来抽取其他项.
energy.xvg和box.xvg文件中给出的是什么性质估计压力,体积和密度的稳定值.
echo"48\n500"|gmxenergy-fener.edr-ocoulomb-inter.xvgecho"49\n510"|gmxenergy-fener.edr-ovanderwaals-inter.xvgcoulomb-inter.xvg和vanderwaals-inter.xvg文件中给出的是什么量
质量保证中最重要的检查事项之一就是确保周期性映像之间没有直接的相互作用.由于周期性映像是全同的,其间的相互作用是物理上不应该发生的自相互作用,会导致模拟结果不合理.设想具有偶极矩的蛋白会有直接的相互作用.那么同一蛋白处于盒子边界的两个末端之间就会产生吸引,这将影响蛋白质的自身的行为并导致模拟结果不合理.我们通过计算每一时刻周期性映像之间的最小距离来证实这种情况不会发生.这通过mindist来实现:
gmxmindist-ftraj.xtc-stopol.tpr-odminimal-periodic-distance.xvg-pi周期性映像之间的最小距离多大,何时出现在你的模拟过程中,用于长程非键作用的截断距离是多少根据这一距离,两个周期性映像之间允许的最近距离多大当最小距离小于截断距离时,哪一非键能量项受影响最大,为什么如果最小距离小于截断距离,会发生什么情况你的模拟中是否出现了这种情况选择C-alpha组重新运行mindist,结果有变化么对你的体系而言,这意味着什么
要注意小距离事件是不时出现还是持续出现.如果持续出现,很可能影响蛋白的动力学;如果只是偶尔出现,那基本没有影响.
不仅直接相互作用需要担心,非直接效应,即以水为媒介的间接相互作用也可能引起问题.例如,蛋白可以导致水在离其最近的四个溶剂化层中出现有序性,这大约对应于1nm的距离.理想情况下,最小距离不应该小于2nm.
gmxgyrate-ftraj.xtc-stopol.tpr-oradius-of-gyration.xvg查看回旋半径及其分量,注意这些值如何达到平衡值.
回旋半径收敛了么如果是,什么时候收敛,收敛值为多少不同多肽模拟给出的回旋半径其涨落行为是否类似,为什么你能将回旋半径的涨落与周期性映像的最小距离联系起来么
gmxrmsf-ftraj.xtc-stopol.tpr-ormsf-per-residue.xvg-oxaverage.pdb-oqbfactors-residue.pdb-res使用xmgrace查看RMSF的图形,区分柔性和刚性区域.
确定最柔性区域的起始和终止残基编号,给出其最大振幅.比较不同多肽构象的结果.是否有区别如果是,哪个构象柔性最大,哪个构象柔性最小
为了对这些结果获得关联性的印象,这里有个人类朊蛋白质1qlz的RMSF,图中标示出了会导致CJD疾病的突变残基.
将两个pdb文件载入PyMOL,根据b因子对结构bfactors.pdb进行着色,并确定柔性区域.平均结构是非物理结构.查看下其中的一些侧链,注意观察平均对构象的影响.
pymolaverage.pdbbfactors-residue.pdbspectrumb,selection=bfactors-residue下图是根据计算的b因子对1KLU多肽的参考结构进行着色后的结果(提示:下图的设置如下:assticks;showspheres;setstick_radius,0.12;setsphere_scale,0.25)
颜色分布在b因子值的范围内,其中蓝色表示最小值(最稳定),红色表示最高值(波动最大).你可以通过截断最大值来调整颜色范围,例如将其设置为350:
Q=350;cmd.alter("all","q=b>QandQorb");spectrumq,selection=bfactors-residue如果你很好奇,可以再次计算b因子,但这次使用每个原子的值:
g_rmsf-ftraj.xtc-stopol.tpr-ormsf-per-atom.xvg-oqbfactors-atom.pdb在当前的PyMOL中载入新生成的bfactors-atom.pdb文件,这样可以直接将其与残基平均的b因子相比较.如果需要,不要忘了重新调整颜色范围.
loadbfactors-atom.pdbspectrumb,selection=bfactors-atom比较和对照两个b因子的结构.
以下图像显示的是根据模拟计算得到人类野生型UbcH8蛋白的b因子着色图.蓝色对应低值,红色对应高值.白色区域表示目前已知的能够反转蛋白质相互作用特征的残基.在图像右侧,你可以看到那些在螺旋2的前后环区有较高的b因子.第二个图像是螺旋2的前环的放大显示.
当心!由于你的多肽可能会跳出盒子外,我们必须处理轨迹,将粒子重新置于中心的周期性映像中.为此,可使用下面的命令:
再次计算,但只考虑骨架原子:
gmxrms-ftraj_nojump.xtc-stopol.tpr-ormsd-backbone-vs-start.xvg这次的RMSD值更低,这是由于计算时排除了通常更柔性的侧链原子,在两种情形下.两个RMSD都应该增大到一个平台值,这意味着相对于参考结构,多肽的结构达到了一定的距离,然后或多或少保持在那个距离.然而,随着距离的增加,可能的构象数目也在增加,这意味着尽管RMSD达到了一个平台值,但结构可能仍在趋近于其平衡态.为此,我们建议同时也检查趋向于平均结构的收敛情况.
echo1|gmxtrjconv-ftraj_nojump.xtc-stopol.tpr-otraj_protein_nojump.xtcgmxrms-ftraj_protein_nojump.xtc-saverage.pdb-ormsd-all-atom-vs-average.xvggmxrms-ftraj_protein_nojump.xtc-saverage.pdb-ormsd-backbone-vs-average.xvg与前面得到的图像进行比较.注意在哪一点RMSD值开始趋平.
简要讨论下两个(相对于起始结构与相对于平均结构)图像的区别,哪一个更能表征收敛情况
到从为止,我们完成了分析的第一部分,包括可视化考察和质量确保检查.现在该深入一点,发掘一下蛋白质内部的情况.分析的第二部分包括可根据多肽构象计算的结构性质,例如氢键数量,溶剂可及表面或特定的原子-原子间距离等.
说明:
氢键数目是一类信息丰富的性质,无论是内部氢键(蛋白-蛋白)还是蛋白和周围溶剂之间的氢键.氢键存在与否可以通过氢键施体-H-受体之间的距离和施体-H-受体之间角度来推断.hbond命令可计算模拟过程中分子间或组间的氢键数目以及氢键距离或角度的分布.使用下面的命令,然后查看得到的输出文件
echo11|gmxhbond-ftraj_nojump.xtc-stopol.tpr-numhydrogen-bonds-intra-protein.xvgecho112|gmxhbond-ftraj_nojump.xtc-stopol.tpr-numhydrogen-bonds-protein-water.xvg讨论两种情形下氢键数目的关系,每种情形下的氢键数目的波动情况.
特定的氢键可以使用包含待研究原子编号的索引文件来考察.查看多肽前一半和后一半所涉及的氢键.你必须看看confout.gro文件来检查残基编号并将多肽大致分为两半.假定你的多肽含有14个残基,其编号始于22终于35.你想将它分为两部分,22-28和29-35.下面的第一个命令将在组选择菜单中创建两个新的条目,part_1和part_2.第二个命令是一个使用索引文件运行hbond的通用命令.你可能想要看看N端和C端之间的氢键,例如,可以用来监测β发卡的形成.
echo"r22-28\nname19part_1\nr29-35\nname20part_2\nq"|gmxmake_ndx-fconfout.gro-omy_index.ndxgmxhbond-ftraj_nojump.xtc-stopol.tpr-numhydrogen-bonds-strand.xvg-nmy_index.ndx注意,你必须根据自己的多肽替换上面命令中的数字.
根据氢键分析,你的多肽在模拟过程中是否形成类似β发卡的构象
GROMACS提供了一个dssp的接口,可以计算轨迹中每帧的二级结构.
首先,需要生成消除跳跃的轨迹
gmxtrjconv-ftraj.xtc-otraj_nojump.xtc-pbcnojump-dt50然后运行以下命令
gmxxpm2ps-fsecondary-structure.xpm-osecondary-structure.eps-size4000ps2pdfsecondary-structure.epsevincesecondary-structure.pdf讨论二级结构的变化,如果有的话.比较不同蛋白二级结构的稳定性.
gmxrama-ftraj_nojump.xtc-stopol.tpr-oramachandran.xvg此文件中包含了全部氨基酸的所有phi/psi角.要抽取某个残基的角度,例如LEU-24,可键入
egrep"@|LEU-24"ramachandran.xvg>rama-LEU-24.xvg抽取每个残基(除边界外)的拉氏图,并用xmgrace对其进行可视化,描述它们的主要相似性与差别.
根据拉氏图,讨论每个残基的构象稳定性.
每组学生模拟了一个从不同构象开始的多肽,这样做的目的在于增加采样,或者增加模拟能够覆盖的结构空间.为了对轨迹和多肽构象性质有一个完整的观点,我们必须将轨迹拼合到一起.
gmxpdb2gmx-fXXXX_peptide.pdb-oXXXX_peptide.gro-ignh这样,使用新生成的这个文件来创建索引文件,这样我们可以进行更多的选择.在这一步,与前面不同,我们不需要指定任何特定的残基组,因此,我们简单地推出程序(键入q)
gmxtrjconv-ftraj_helix.xtc-otraj_helix_10-50ns.xtc-nXXXX_peptide_index.ndx-pbcnojump-dt100-b10000为什么我们只分析最后40ns的轨迹
gmxtrjcat-ftraj_bound_10-50ns.xtctraj_helix_10-50ns.xtctraj_polypro_10-50ns.xtctraj_extended_10-50ns.xtc-otraj_concat_10-50ns.xtc-cat-settime5.又是RMSD我们在前面已经计算过RMSD,并用其来检查模拟的收敛性,但它也可以用于更进一步的分析.RMSD是两个结构之间的表征.如果我们对轨迹文件中每一对结构的组合计算RMSD,就可以看到是否有属于同一类型或具有相同特征的结构组.属于同一组的结构其RMSD值较低,而与其他组结构的RMSD值更高.利用矩阵来表示RMSD值,可以用于识别转变状态.
要建立RMSD矩阵,可使用rms处理两条轨迹.如果你要单独考虑组(簇)及其在不同轨迹之间的转变,可以将所有的四条轨迹合并为一条,再使用rms生成交叉RMSD矩阵.所有的RMSD计算时,选择主链组Mainchain.
gmxrms-sXXXX_peptide.gro-ftraj_concat_10-50ns.xtc-f2traj_concat_10-50ns.xtc-mrmsd-matrix.xpm-tuns得到的矩阵是灰阶图.为了显示更清楚,可以使用彩虹梯度图.
gmxxpm2ps-frmsd-matrix.xpm-ormsd-matrix.eps-rainbowblueps2pdfrmsd-matrix.epsevincermsd-matrix.pdf分析时为什么选择主链组你看到多少簇在不同的轨迹中,你是否对相同构象进行了采样你是如何发现的什么是轨迹的重叠,即对相同构象空间进行采样你能发现多少转变从这个分析你能得到什么结论这是你预期的结果么请证实你的观点.
基于结构间的距离RMSD,可以将结构归并为反映构象可及性范围及其相对权重的一组组团簇.这可以通过聚类算法来完成,cluster命令实现了一些聚类算法.这一程序会生成好几个输出文件.检查该程序的帮助文档,了解每一个它们的含义,然后运行程序.注意,我们已经计算了RMSD矩阵,可以将它作为cluster的输入.
gmxcluster-hgmxcluster-sXXXX_peptide.gro-ftraj_concat_10-50ns.xtc-dmrmsd-matrix.xpm-distrmsd-distribution.xvg-oclusters.xpm-szcluster-sizes.xvg-trcluster-transitions.xpm-ntrcluster-transitions.xvg-clidcluster-id-over-time.xvg-clclusters.pdb-cutoff0.2-methodgromos聚类算法使用的RMSD截断值为多大这一值代表什么,如何影响得到的团簇数目共有多少团簇最大的两个其尺寸多大
更改RMSD的截断值,以获得适当数目的团簇(10到100之间).使用PyMOL打开cluater.pdb,比较前两个团簇的结构.
disableallsplit_statesclustersdeleteclusters/foriinrange(3,100):cmd.delete("clusters_%04d"%i)dssshowcartoonutil.cbamclusters_0002alignclusters_0001andssh,clusters_0002andssh前两个团簇之间是否存在可觉察的差别
采用RMSD比较结构的一个不足之处在于,它包含了最小二乘叠合,这会影响结果.但是,一个蛋白的结构也可以使用一系列的原子间距来表示.这可以用于获得一个比较的表征量,并且不依赖于叠合,这就是距离RMSD(dRMSD).可以利用rmsdist命令来计算dRMSD.
gmxrmsdist-sXXXX_peptide.gro-ftraj_concat_10-50ns.xtc-odistance-rmsd-all-atoms.xvgdRMSD何时收敛收敛值多大所得图像与标准RMSD相比有何区别查阅GROMACS手册,看二者是如何计算的.试着根据二者的计算方法解释你观察到的区别.
我们现在要回答构象采样的问题,在我们的轨迹中,我们是否对多肽的结合构象进行了采样为此,我们要对多肽的结合构象再次运行rms.我们可以使用合并的轨迹来进行这一分析,而不是计算并生成四个不同的RMSD图形.记住你生成的合并轨迹中构象的顺序,这很重要.当叠合计算RMSD值时,选择主链Mainchain组.
gmxrms-ftraj_concat_10-50ns.xtc-sXXXX_peptide.gro-ormsd_concat-mainchain-vs-bound.xvgxmgracermsd_concat-mainchain-vs-bound.xvggmxrmsdist-ftraj_bound_10-50ns.xtc-sXXXX_peptide.gro-odist-rmsd_bound-mainchain-vs-bound.xvgxmgracedist-rmsd_bound-mainchain-vs-*.xvg如果在计算时我们包含整个蛋白(即所有原子),而不是只选择主链Mainchain原子,结果会有什么变化RMSD和距离RMSD图形有区别么在不同轨迹中,你是否对多肽的结合构象进行了采样对你研究的情况,关于对蛋白多肽识别的构型采样假说你能得出什么结论
本分析的最后一步是计算自由能形貌(FEL,freeenergylandscape),它由不同多肽构象开始的轨迹采样所得.除了看起来很酷以外,它能显示出在整个模拟过程中多肽所经历的自由能的谷和丘.此外,在对接计算中需要选择有代表性的结构时,FEL会变得更有意义.
FEL表示一个映射,分子在模拟过程中所经历的所有可能构象到相应能量的映射.典型的能量是Gibbs自由能.正如你所想象的,使用直角坐标来代表不同的构象是不合适的(你能想象4维形貌么).因此,FEL通常使用两个变量来表示,它们反映了体系的特定性质,并表征了构象变化.例如你可以使用绕一根特定键的扭转角,或分子的回旋半径,或相对于天然状态的RMSD来作为这两个变量.第三个变量是自由能,可以从体系相对前面所选变量的分布(概率分布)来估计.当使用三维表示时,形貌图中的谷表示低自由能区域,代表体系的亚稳定构象,丘表示连接这些亚稳定状态的能量势垒.
我们如何从MD轨迹中抽取代表性结构用于对接模拟给一个前面用来分析多肽构象变化的方法,并说明这种方法可用于此目的.
因为FEL分析依赖于在多肽构象空间的充分采样,并基于构象的分布进行Gibbs自由能估计,我们需要准备新的轨迹文件(.xtc),其中的帧数是我们到目前为止所用的十倍.注意我们没有使用-dt选项.对所有四条原始轨迹重复下面的步骤,注意要更改名字.
gmxtrjcat-ftraj_bound_10-50ns-highresolution.xtctraj_helix_10-50ns-highresolution.xtctraj_polypro_10-50ns-highresolution.xtctraj_extended_10-50ns-highresolution.xtc-otraj_concat_10-50ns-highresolution.xtc-cat-settime现在我们需要生成FEL计算所需的数据.先计算合并轨迹的回旋半径Rg,方法如前:
gmxgyrate-ftraj_concat_10-50ns-highresolution.xtc-sXXXX_peptide.gro-org-concatenated_traj.xvg为计算相对于平均结构的RMSD,首先必须重新计算整条轨迹的平均结构.再次使用rmsf命令生成平均结构(-ox选项),然后使用平均结构做参考,利用rms计算RMSD,方法如前.当提示时,对叠合和计算都选择蛋白Protein组.
perlsham.pl-i1rg-concatenated_traj.xvg-i2rmsd-vs-average-concat.xvg-data11-data21-ogsham_input.xvg现在我们有了正确的sham输入文件,可用于生成FEL.如果你还记得,Gibbs自由能可以根据分布概率计算出来,并且依赖于指定的温度.使用sham的选项-tsham来指定正确的体系温度(如果你不记得这个值,可检查成品模拟所用mdp文件中的值).
gmxsham-fgsham_input.xvg-lsfree-energy-landscape.xpm如果使用sham计算FEL时,指定非常高的温度,会出现什么情况
从FEL中选择你认为能代表多肽状态的五个点,记下它们的坐标.
pythonget_timestamp.py-fgsham_input.xvg-10.657-23.1415如果你总是得到Notimestampfound...这样的错误信息,试着使用附近稍有区别的点,或打开脚本增加变量nval的值到100或更大.
gmxtrjconv-ftraj_concat_10-50ns-highresolution.xtc-sXXXX_peptide.gro-dump45000-orepresentative_45000ps.pdb在PyMOL中打开代表性结构以及结合结构,并比较二者.
以下分析整理自泛素耦合酶教程,用于说明操作过程,但具体问题未必适用于前面的多肽-蛋白教程.
特定氢键可以用包含相应原子数量的索引文件来得出.从分析1得出的RMSF及b-factors显示loops2和3(helix2附近)值比较高.实验数据也显示,loop1可能在UbcH6andUbcH8的多个行为中起了一定作用.看看这些loop所包含的氢键连接.第一个命令会在菜单中弹出3个新的条目去选择基团,每个loop一个.第二个命令是一个通用命令g_hbond,需要一个索引文件.你可能想看看每个loop里的氢键、loop1和loop2之间的氢键；例如,某个特定的loop和蛋白质其他部分的氢键(为此可能需要修改索引文件)或者和水形成的氢键,等…
mkdirsaltbridgecdsaltbridgegmxsaltbr-f../traj_nojump.xtc-s../topol.tpr-t0.5-sep为了更清晰,删除与钠离子和氯离子有关的文件:
rm-f*CL-**NA+*\#*看看以下残基之间的相互作用:
对于这些相互作用你有什么想法残基K60和D88高度保守,仅仅将UbcH8中的残基D56突变为E56都会改变蛋白的相互作用方式.其可能原因何在
gmxsasa-ftraj_nojump.xtc-stopol.tpr-osolvent-accessible-surface.xvg-oaatomic-sas.xvg-orresidue-sas.xvg哪个残基是最容易被溶剂触及的
前面用到的rmsdist命令也可用于进一步的距离分析.特别地,为了解结构及其稳定性,查看原子间的平均距离及其波动可能会有用处.使用下面的命令对蛋白质中每对原子间的平均距离及其波动进行计算获得矩阵,然后用上面的步骤重新着色,并显示.png文件的结果(rmsmean.png和rmsdist.png).
gmxrmsdist-stopol.tpr-ftraj.xtc-rmsrmsdist.xpm-meanrmsmean.xpm-dt10简单地解释两个图像:rmsmean表示结构,rmsdist表示柔性/稳定性.回想前面分析得到的信息并查看结构.
gmxrmsdist-stopol.tpr-ftraj.xtc-nmr3nmr3.xpm-nmr6nmr6.xpm-noenoe.dat-dt10给nmr3.xpm和nmr6.xpm重新着色,查看这些矩阵,也查看noe.dat文件.
模拟战,最小的1/r3和1/r6平均距离为多少
N-H序参量可以用chi命令计算.这个程序可以写出一个.pdb文件,把序参量加入了b因子列,更容易查看.该程序也计算J耦合参数,并可以和NMR结果相比较,或用于指导NMR实验.
gmxchi-stopol.tpr-ftraj.xtc-oorder-parameters.xvg-porder-parameters.pdb-jcJcoupling.xvg看看.xvg文件中的序参量,并用PyMOL查看.pdb文件,根据b因子的值给残基着色.
记下起始与终止的残基,具有最高序参量值的两个区域的平均值.序参量与波动(RMSF)相比如何
一个常用的,但常常理解不深的分析方法是对轨迹进行主成分分析(PCA,PrincipalComponentsAnalysis).这种方法有时候也称为本征动力学(ED,essentialdynamics),目的在于识别原子的大尺度集约运动,从而揭示隐藏在原子波动后面的结构信息,帮助确定哪些运动方式对蛋白质的整体动力学贡献最大.
协方差分析会生成很多文件,因此最好在一个新的目录中运行:
mkdirCOVARcdCOVAR协方差矩阵的构建和对角化可使用covar命令.键入下述命令进行分析:
协方差矩阵的维数多少特征值的总和多少
现在,看看covar.xpm文件中的协方差矩阵.
xviewcovar.xpm矩阵显示了原子间的协方差.红色表示两个原子运动一致,而蓝色表示它们彼此向相反的方向运动.红色的深度表示波动振幅的大小.对角线上的值对应于前面得到的RMSF图.
查看除端基外运动最剧烈的两个部分,它们之间如何相对运动,它们各自相对于蛋白的其他部分如何运动
使用文本编辑器查看eigenvalues.xvg文件,计算头五个特征向量对总运动的百分比以及累积百分比.
典型地,最初五个特征值将捕获主要运动,这相当于>80%的总运动.如果解释的总波动较低,就说明没有明确的集约运动.
为了解特征值的实际意义,可用anaeig命令作进一步的分析.为了更近地看看前两个特征值,键入以下命令
gmxanaeig-s../topol.tpr-f../traj.xtc-veigenvectors.trr-eigeigenvalues.xvg-projproj-ev1.xvg-extrev1.pdb-rmsfrmsf-ev1.xvg-first1-last1gmxanaeig-s../topol.tpr-f../traj.xtc-veigenvectors.trr-eigeigenvalues.xvg-projproj-ev2.xvg-extrev2.pdb-rmsfrmsf-ev2.xvg-first2-last2特征值对应于运动方向.选项-extr沿着选定的特征值从轨迹中提取极端结构.把这些结构导入PyMOL查看:
pymolev.pdb把pdf文件中的模型分开,删除原始结构.
split_statesev1split_statesev2deleteev1orev2给模型着色.特征值1中的极端结构显示为蓝-绿色而特征值2为黄-红色.
spectrumcountdsshideeverythingshowcartoon用PyMOL的align命令,能画出表示两种构象差异的小条.
alignev1_0001and(n.c,n,ca),ev1_0002and(n.c,n,ca),object=diff1alignev2_0001and(n.c,n,ca),ev2_0002and(n.c,n,ca),object=diff2对特征值1而言,极端结构之间的最大区别是什么对特征值2呢
为了理解特征值的意义,想象一下旅行推销员在欧洲城市间的移动.这种移动可用地球坐标系统来说明,对每个位置需要采用三个坐标.虽然这样做没有问题,但如果你只想解释推销员的移动,这种方法不是最佳的.因为理论上,任何坐标系统都一样好,我们可以定义一个新的坐标系统来解释推销员的移动.实际上,因为地球表面可以用二维空间代表,我们只需要两个坐标而无须三个.直观看来,有人会以南北极(经度)和东西轴(维度)说明,但也可以从移动中推断出轴.比方说推销员在伦敦-雅典轴上走的最多,这个轴可以作为第一个特征值;第二个特征值与第一个正交.用这种方式,推销员在欧洲的每个位置就可以用在这两个特征值上的投影唯一地确定下来.对它们的投影作图,就可以显示出旅行路线.沿着第一个坐标的极端投影对应于雷克雅未克(Reykjavik,冰岛首都)和莫斯科,即使它们实际上不在这个轴上.
蛋白质构象也和上面所述的一样.你看到的极端投影并不一定对应于物理结构,但是它们可以表征沿轴的运动和总的运动程度.为了解蛋白质沿构象空间的移动,我们可以画出特征值2对特征值1的投影图.为此,从两个.xvg文件中提取投影数据并且合并到文件ev1-vs-ev2.dat中.注意’>’表示输出由屏幕重定向到一个文件中,所以你看不到任何屏幕输出.
grep-v"^[#@]"proj-ev1.xvg|awk'{print$2}'>proj-ev11.datgrep-v"^[#@]"proj-ev2.xvg|awk'{print$2}'>proj-ev12.datpasteproj-ev11.datproj-ev12.dat>ev1-vs-ev2.dat为了解释模拟过程,我们也提取最后7.5ns(最后1500点)和5.0ns(最后1000点)的投影,然后将它们导入xmgrace.
tail-n1500ev1-vs-ev2.dat>last-7.5ns.dattail-n1000ev1-vs-ev2.dat>last-5.0ns.datxmgraceev1-vs-ev2.datlast-7.5ns.datlast-5.0ns.dat投影的形状如何它们相互依赖么(分布有所重叠)如果只对最后7.5ns进行分析,是否得到相同的特征向量(轴)使用最后5.0ns呢
Writeaconcludingparagraph,comparingtheresultsobtainedforthedifferentproteins.AlsoreflectontheoverallstabilityandtheprobabilitythatthestructuredepositedinthePDBproperlyreflectsthesolutionstructure.Inotherwords,doesthestructurestayclosetothestartingstructureordoesitdriftawayandhowmuchCanyousuggestamechanismthatwouldexplainthedifferencesinUbcH6andUbcH8interactionprofilesandwhyasingleconservedmutationcanrestoretherichinteractionprofileofUbcH6whenstartingfromUbcH8Whatarethelimitsofthesimulationstoexplainsuchacomplexbiologicalprocess